Pointwise entropy in multiple alignments

Shannon entropy is used to quantify the variability of positions in a multiple sequence alignment. This application is particularly widespread due to the ubiquitous use of sequence logos [15, 16] to present motifs in proteins, DNA and RNA. Letting denote the 4-letter alphabet {A, C, G, U}, the pointwise entropy H₁(k) at position k in the alignment is defined by , where p_a is the proportion of nucleotide a at position k. Entropy values range from 0 to log 4, where high entropy entails uncertainty or disagreement of the nucleotides at position k. Average pointwise sequence entropy is often expressed in bits, where logarithm base 2 is used instead of the natural logarithm. The concept of sequence logo has many generalizations; indeed, logos for DNA major groove binding are described in [16], logos for tertiary structure alignment of proteins are described in [17], logos for RNA alignments including mutual information on base pair covariation are described in [18], and logos with secondary structure context of RNAs that bind to specific riboproteins are described in [19, 20].

Positional entropy

For a given RNA sequence , and for 1 ≤ i < j ≤ n, define the base pairing probability p_i,j to be the sum of Boltzmann factors of all secondary structures that contain base pair (i, j), divided by the partition function, i.e. (3)

Here p(s) is the Boltzmann probability of structure s of , E(s) is the Turner free energy of secondary structure s[5], R ≈ 0.001987 kcal/mol.K is the universal gas constant, T is absolute temperature, and the partition function Z = ∑_s exp(−E(s)/RT), where the sum is taken over all secondary structures s of a. Base pairing probabilities can be computed in cubic time by McCaskill’s algorithm [11], as implemented in various software, including the Vienna RNA Package RNAfold -p[10].

Define the positional base pairing probability distribution at fixed position 1 ≤ i ≤ n by (4)

For each fixed value of i, is a probability distribution, where j ranges over 1,…,n, the positional structural entropy H₂(i) at position i is defined by (5)

Low values of positional entropy at position i indicate that there is a strong agreement among low energy structures in the Boltzmann ensemble that either i is unpaired, or that i is paired with the same position j. The average positional entropy 〈H₂〉 is the average taken over all positions of the sequence. Structural positional entropy was first defined by Huynen et al. [21], who used the term S-value for average positional entropy, and showed that RNA nucleotide positions having low entropy correspond to positions where the minimum free energy (MFE) structure tends to agree with that determined by comparative sequence analysis. In [22], Mathews made a similar analysis, where in place of S-value, a normalized pseudo-entropy value was used, defined by −∑_{1 ≤ i < j ≤ n} p_i,j ln p_i,j/n. Positional entropy of RNA secondary structures can be presented by color-coding each nucleotide, where the color of the kth nucleotide reflects the positional entropy H₂(k) as defined in Eq (5). The Rfam 12.0 database [23] uses such color-coded secondary structures, since the base-pairing of positions having low entropy is likely to be correct [21, 22].

Derivational entropy using stochastic context free grammars

Manzourolajdad et al. [24], Sukosd et al. [25] and Anderson et al. [26] describe the computation of structural entropy for stochastic context free grammars (SCFGs), defined by −∑_s p(s) ln p(s), where the sum is taken over all secondary structures s of a given RNA sequence, and p(s) is the probability of deriving the structure s in a particular grammar G, defined as follows. Suppose that S = S₀ is the starting nonterminal for the grammar G, s = S_m is the secondary structure s consisting only of terminal symbols belonging to the alphabet , and that S₁,…,S_m−1 are expressions consisting of a mix of nonterminal and terminal symbols. If S₀ →_G S₁ →_G S₂ →_G⋯ →_G S_m is a leftmost derivation using production rules from grammar G and for each i = 0, …, m−1, we let p_i denote the probability of applying the rule S_i → S_i+1, then p(s) is defined to be the product . It should be noted that the derivational probability p(s) heavily depends on the choice of grammar G as well as on the rule application probabilities p_i, obtained by applying expectation maximization to a chosen training set of secondary structures.

Anderson et al. [26] are motivated to compute derivational entropy of a multiple alignment of RNAs, in order to provide a numerical quantification for the quality of the alignment—specifically, their paper shows that accurate alignment quality corresponds to low derivational entropy. In [27], Sukosd et al. describe the software PPfold, a multithreaded version of the Pfold RNA secondary structure prediction algorithm. Subsequently, Sukosd et al. [25] describe how to compute the derivational entropy for the grammar used in the PFold algorithm (grammar G6 as defined in [28]), and show that derivational entropy is correlated with the accuracy of PPfold structure predictions, as measured by F-scores. In contrast, Manzourolajdad et al. [24] computed the derivational entropy of various families of noncoding RNAs, using the trained stochastic context free grammars G4,G5,G6 [28], which they denote respectively as RUN (G4), IVO (G5) and BJK (G6). The Linux executable and trained models can be downloaded from http://rna-informatics.uga.edu/malmberg/ for three RNA stochastic context free grammars, each with three trained models using the training sets ‘Rfam5’, ‘Mixed80’, and ‘Benchmark’—see [24] for description.

The plan of the remainder of this paper is as follows. In the Methods, we provide a description of our novel entropy algorithms, beginning with an overview in Section “Statistical mechanics”, where we derive a relation between structural entropy H and expected energy 〈E〉. This relation allows us to provide a crude estimate of H by sampling. Expected energy can be computed from the derivative of the logarithm of the partition function with respect to temperature; a finite difference computation then yields our first algorithm to compute structural entropy, while a dynamic programming approach for the expected energy yields our second algorithm. In the Results, we compare our structural entropy software, RNAentropy, with software for SCFG derivational entropy, and then use RNAentropy in several applications. Section “Comparison of structural entropy and derivational entropy” benchmarks the time required to compute structural entropy using our two algorithms, versus the time required to compute derivational entropy using the program of [24]. Numerical values for structural and derivational entropies are compared, along with their distributions. In addition to a comparison of run times and entropy values, we compute the correlation of structural entropy, derivational entropy, and a variety of measures, such as ensemble defect [29], positional entropy [21], structural diversity [30], etc. Such measures have recently been used in the design of experimentally validated RNA molecules [14, 31]. Motivated by the fact that calmodulin-ligand binding has been shown to depend on conformational entropy [1], in Section “Correlation with hammerhead cleavage activity”, we show an improvement in the correlation between hammerhead ribozyme cleavage activity and total change of energy [32], if conformational entropy is also taken into account. In Section “Structural entropy of HIV-1 genomic regions”, we compute the entropy of genomic portions of the HIV-1 genome and compare entropy Z-scores with known HIV-1 noncoding elements. Finally, in the Discussion, we describe differences between the methods and discuss the numerical discrepancy between thermodynamic structural entropy values and SCFG derivational entropy values. For more background on RNA, an excellent, though somewhat outdated, review of computational and physical aspects of RNA is given by Higgs [33].

Statistical mechanics

Shannon entropy for the Boltzmann ensemble of secondary structures of a given RNA sequence is defined by (6) where G denotes the ensemble free energy −RT ln Z. It follows that if the energy E(s) of every structure s is zero, or if the temperature T is infinite, then entropy is equal to the logarithm of the number of structures. Note as well that in the Nussinov energy model [34], where each base pair has an energy of −1, it follows that the expected energy is equal to −1 times the expected number of base pairs, i.e. 〈E〉 = −∑_{i < j} p_i,j, where p_i,j is the probability of base pair (i, j) in the Nussinov model.

By sampling RNA structures with the RNAsubopt program from Vienna RNA Package [10], we can approximate the value of expected energy, and hence obtain an approximation of the thermodynamic entropy by using Eq (6). This can be done in two distinct manners.

In the first approach, a user-specified number N of low energy structures from the thermodynamic ensemble can be sampled by using the algorithm of Ding and Lawrence [35], as implemented in RNAsubopt -p N. A sampling approximation for the expected energy is then defined to be the arithmetic average of the free energy of the N sampled structures. In the second approach, all structures can be generated, whose free energy lies within a user-specified range E of the minimum free energy, by using the algorithm of Wuchty [36], as implemented in RNAsubopt -e E. Let Z₀ be an approximation of the partition function, defined by summing the Boltzmann factors exp(−E(s)/RT) for all generated structures. Define the (approximate) Boltzmann probability of a generated structure s to be p(s) = exp(−E(s)/RT)/Z₀. An approximation for the expected energy is in this case taken to be ∑_s p(s) ⋅ E(s), where the sum is taken over all structures s, whose free energy is within E kcal/mol of the minimum free energy. In either case, the resulting entropy approximation is not particularly good. For instance, the thermodynamic entropy of the 78 nt arginyl-tRNA from Aeropyrum pernix (accession code tdbR00000589 in the Transfer RNA database tRNAdb [37]) is 5.44, as computed by the algorithm RNAentropy described in this paper, while the entropy approximation by the first sampling approach with N = 10,000 is 4.71 and that of the second sampling approach with E = 10 is 4.68. Since the estimate from each sampling approach has greater than 13% relative error, sampling cannot be used to provide accurate entropy values. For that reason, we now briefly describe two novel, cubic time algorithms to compute the exact value of structural entropy—details of the algorithms are further described in Sections “Entropy by statistical physics” and “Entropy by dynamic programming”.

Algorithm 1: Formal temperature derivative (FTD)

It is well-known from statistical physics that the average energy 〈E〉 of N independent and distinguishable particles is given by the following formula (cf equation (10.36) of [38]): (7)

This equation does not hold in the case of RNA secondary structures with the Turner energy model; however, Eq (7) is close to being correct. The idea of Algorithm 1 is to use finite differences to approximate the derivative , thus obtaining the expected energy 〈E〉, from which we obtain the structural entropy by applying Eq (6). As shown later, certain technically subtle issues arise in this approach; in particular, the derivative must be taken with respect to the formal temperature, which represents only those occurrences of the temperature variable within the expression RT. Formal temperature is distinct from table temperature, which latter designates all occurrences of the temperature variable in the Turner energy parameters. This will be fully explained in Section “Entropy by statistical physics”. For this reason, Algorithm 1 is named FTD, for formal temperature derivative.

Algorithm 2: Dynamic Programming (DP)

Recall that the partition function for a given RNA sequence a is defined by Z = ∑_s exp(−E(s)/RT), where the sum is taken over all secondary structures of a. Letting BF(s) = exp(−E(s)/RT) denote the Boltzmann factor of s, it follows that the Boltzmann probality of secondary structure s satisfies p(s) = BF(s)/Z, and hence (8) where Q = ∑_s BF(s) ⋅ E(s). The partition function Z can be computed by McCaskill’s algorithm [11], while in Sections “Entropy by dynamic programming”, we describe a dynamic programming algorithm to compute Q(a). Since this method uses dynamic programming, Algorithm 2 is named DP.

Both FTD and DP support the Turner’99 and Turner’04 energy models [5], and all references to FTD and DP mean FTD’04 and DP’04, unless otherwise stated (there are small numerical differences in the entropy, depending on the choice of Turner parameters). Moreover, both algorithms allow the user to specify an arbitrary temperature T for the computation of structural entropy. This latter feature could prove useful in the investigation of thermoswitches, also called RNA thermometers, discussed later. The software RNAentropy implements both algorithms, and is available at http://bioinformatics.bc.edu/clotelab/RNAentropy.

Entropy by statistical physics

Here we show that for the Turner energy model of RNA secondary structures, expected energy satisfies (9) although equality does not strictly hold. Indeed, (10)

Let formal temperature denote each occurrence of the temperature variable T within the expression RT, while table temperature denotes all other occurrences (i.e. table temperature refers to the temperature-dependent Turner free energy parameters [5]). This will shortly be explained in greater detail. From Eq (10), it follows that expected energy 〈E〉 is equal to RT² times the derivative of ln Z(T) with respect to formal temperature, which latter we define to be the formal temperature derivative of ln Z(T).

If we treat the energy E(s) of structure s as a constant (computed at either the default temperature of 37°C, or at a user-specified temperature T), then the second term of Eq (12) disappears, and we can approximate by the finite difference , where for instance ΔT = 10⁻⁷. This requires a modification of McCaskill’s algorithm [11] for the partition function Z(T), where we distinguish between formal temperature and table temperature. Our software RNAentropy implements such a modification, and thus supports the formal temperature derivative (FTD) method of computing thermodynamic structural entropy.

Note that the function ln Z(T) is decreasing and concave down, so barring numerical precision errors, the finite difference is negative and slightly larger in absolute value than the formal temperature derivative . From Eq (6), structural entropy H is equal to 〈E〉/RT + ln Z and so there will be a small numerical deviation between the value of H, computed by the FTD (formal temperature derivative) method currently described, and the exact value of H computed by the DP (dynamic programming) method, described in Section “Entropy by dynamic programming”. In particular, entropy values computed by FTD should be slightly smaller than those computed by DP, where the discrepancy will be visible only for large sequence length. This is indeed observed in Fig 1B and in data not shown.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 1.

(A) Average run times, with (tiny) error-bars of ±1 standard deviation, for each of the five methods DP, FTD (ΔT = 10⁻⁷), SCFG(G6,Benchmark), SCFG(G4,Benchmark), and SCFG(G5,Benchmark). Averages were determined for 100 random RNA sequences of length n, each having expected compositional frequency of 0.25 for A,C,G,U, where n ranges from 20 to 500 with increments of 5. Methods tested are as follows: (1) DP: dynamic programming computation of expected energy 〈E〉 and partition function to yield H = 〈E〉/RT + ln Z, with Turner 2004 energy parameters. (2) FTD: formal temperature derivative method which computes , where the temperature increment T+ΔT is applied only to occurrences of T within the expression RT—i.e. formal temperature, as explained in the text. Increment ΔT is 10⁻⁷, and Turner 2004 energy parameters are used. (3) SCFG: computation of derivational entropy using the method of [24], for the grammars G4, G5, G6 with grammar rule probabilities from ‘Benchmark’ data (see [24, 28]). SCFG executables and models downloaded from http://rna-informatics.uga.edu/malmberg/. The methods, ordered from fastest to slowest, are as follows: FTD, DP, G5, G6, G4, where FTD and DP are approximately equally fast, while the slowest methods, G6 and G4, have almost identical run times. DP and FTD are an order of magnitude faster than G6. (B) Average entropy values, with error bars of ±1 standard deviation, computed by the methods DP, FTD (ΔT = 10⁻⁷), SCFG(G4,Benchmark), SCFG(G5,Benchmark), and SCFG(G6,Benchmark) for the same data set as in the left panel. The methods, ordered from those returning smallest entropy values to largest, are as follows: FTD, DP, G6, G4, G5. FTD and DP return essentially identical values, with a small deviation for larger sequences due to the finite approximation of the formal temperature derivative.

https://doi.org/10.1371/journal.pone.0137859.g001

We now show that the expression, , occurring as the second term in the last line of Eq (10), is equal to −T ⋅ 〈S_t〉 where 〈S_t〉 denotes the expected change in entropy using the Turner parameters [5], determined as follows. From statistical physics, the free energy E(s) of a secondary structure s satisfies (11) where H_t(s) [resp. S_t(s)] denotes change in enthalpy [resp. entropy] from the empty structure to structure s using the Turner parameters. The term S_t measures the entropic loss due to stacked base pairs, hairpins, bulges, internal loops and multiloops using parameters obtained from least-squares fitting of UV absorption data. In the Turner energy model, entropy S_t and enthalpy H_t are assumed to be independent of temperature, so it follows from Eq (11) that , and hence (12)

To compute S_t(s) for a given secondary structure s of an RNA sequence a, determine the the free energy E(s, 37) [resp. E(s, 38)] of structure s at 37°C [resp. 38°C] by using Vienna RNA Package RNAeval[10]; it then follows from Eq (11) that S_t(s) = E(s, 37)−E(s, 38). Throughout this paragraph, the reader should not confuse the notion of conformational entropy from Eq (1), which is always non-zero and is computed by the novel algorithms described in this paper, with the notion of Turner change of entropy S_t(s) of secondary structure s, which is always negative due to entropic loss in going from the empty structure to a fixed structure s. Nor should the reader confuse the notion of structural entropy, denoted by H and defined in Eq (2), with Turner change of enthalpy H_t(s) of secondary structure s.

Entropy by dynamic programming

Throughout this section, denotes an arbitrary but fixed RNA sequence. Below, we give recursions for Q(a), defined by Q(a) = ∑_s BF(s) ⋅ E(s), where the sum is taken over all secondary structures s of RNA sequence a, E(s) is the free energy of s, using the Turner 2004 parameters, BF(s) = exp(−E(s)/RT) is the Boltzmann factor of structure s, where R is the universal gas constant and T the temperature in Kelvin.

Recursions are also given for the partition function Z(a) = ∑_sexp(−E(s)/RT), where the sum is taken over all secondary structures of a. It follows that the expected energy (13)

For 1 ≤ i ≤ j ≤ n, the collection of all secondary structures of is denoted . In contrast, if s is a secondary structure of a₁, …, a_n, then s[i, j] is the restriction of s to the interval [i, j], defined by s[i, j] = {(x, y): i ≤ x ≤ y ≤ j, (x, y) ∈ s.

Initial steps

For notational convenience, we define Q_i,i−1 = 0 and Z_i,i−1 = 1. If i ≤ j < i + 4, then for any secondary structure s, the restriction s[i,j] is the empty structure, denoted by j − i + 1 dots with zero energy, and so Q_i,j = 0. As well, the only secondary structure on [i,j] is the empty structure, so Z_i,j = 1.

Now assume that i + 4 ≤ j. Since (14) we treat each sum in a separate case. Let bp(k, j) be a boolean valued function with the value 1 if k can base-pair with j; i.e. a_k a_j ∈ {AU, UA, CG, GC, GU, UG}. For secondary structure , let bp(k, j, s) be a boolean function with value 1 if it is possible to add the base pair (k, j) to s and obtain a valid secondary structure; i.e. without creating a base triple or pseudoknot.

Case 1: j is unpaired in [i, j]. For in which j is unpaired, s = s[i, j − 1], BF(s) = BF(s[i, j − 1]), and E(s) = E(s[i, j−1]). The contribution to Q_i,j in this case is given by Q_i,j−1.

Case 2: j is paired in [i, j]. The contribution to Q_i,j in this case is given by (15)

Putting together the contributions from both cases, we have (16)

Recursions for the Turner nearest neighbor energy model

In the nearest neighbor energy model [5, 39], free energies are defined not for base pairs, but rather for loops in the loop decomposition of a secondary structure. In particular, there are stabilizing, negative free energies for stacked base pairs and destabilizing, positive free energies for hairpins, bulges, internal loops, and multiloops.

In this section, free energy parameters for base stacking and loops are from the Turner 2004 energy model [5]. As in the previous subsection, Q, Z are defined, but now with respect to the Turner model.

(17)

It follows that Z = Z_1,n is the partition function for secondary structures (the Boltzmann weighted counting of all structures of ) and (18)

To complete the derivation of recursions, we must define QB_i,j and ZB_i,j for the Turner model.

To provide a self-contained treatment, we recall McCaskill’s algorithm [11], which efficiently computes the partition function. For RNA nucleotide sequence , let H(i,j) denote the free energy of a hairpin closed by base pair (i, j), while IL(i, j, i′, j′) denotes the free energy of an internal loop enclosed by the base pairs (i, j) and (i′, j′), where i < i′ < j′ < j. Internal loops comprise the cases of stacked base pairs, left/right bulges and proper internal loops. The free energy for a multiloop containing N_b base pairs and N_u unpaired bases is given by the affine approximation a + bN_b + cN_u.

Definition 1 (Partition function Z and related function Q)

• Z_i,j = ∑_sexp(−E(s)/RT) where the sum is taken over all structures .
• ZB_i,j = ∑_sexp(−E(s)/RT) where the sum is taken over all structures which contain the base pair (i,j).
• ZM_i,j = ∑_sexp(−E(s)/RT) where the sum is taken over all structures which are contained within an enclosing multiloop having at least one component.
• ZM1_i,j = ∑_sexp(−E(s)/RT) where the sum is taken over all structures s ∈ Q[i,j] which are contained within an enclosing multiloop having exactly one component. Moreover, it is required that (i,r) is a base pair of x, for some i < r ≤ j.
• Q_i,j = ∑_s E(s)⋅exp(−E(s)/RT) where the sum is taken over all structures s ∈ ss[i,j].
• QB_i,j = ∑_s E(s)⋅exp(−E(s)/RT) where the sum is taken over all structures s ∈ ss[i,j] which contain the base pair (i,j).
• QM_i,j = ∑_s E(s)⋅exp(−E(s)/RT) where the sum is taken over all structures s ∈ ss[i,j] which are contained within an enclosing multiloop having at least one component.
• QM1_i,j = ∑_s E(s)⋅exp(−E(s)/RT) where the sum is taken over all structures s ∈ ss[i,j] which are contained within an enclosing multiloop having exactly one component. Moreover, it is required that (i,r) is a base pair of s, for some i < r ≤ j.

For j − i ∈ {0, 1, 2, 3}, Z(i,j) = 1, since the empty structure is the only possible secondary structure. For j−i > θ = 3, we have (19) (20) (21) (22)

Base Case: For j − i ∈ {−1, 0, 1, 2, 3}, Q_i,j = QB_i,j = 0, Z_i,j = 1, ZB_i,j = ZM_i,j = ZM1_i,j = 0.

Inductive Case: Assume that j − i > 3.

Case A: (i,j) closes a hairpin.

In this case, the contribution to QB_i,j is given by (23)

Case B: (i,j) closes a stacked base pair, bulge or internal loop, whose other closing base pair is (ℓ,r), where i < ℓ < r < j.

In this case, the contribution to QB_i,j is given by the following (24) (25)

In the summation notation , if upper bound b is smaller than lower bound a, then we intend a loop of the form: FOR i = b downto a.

Case C: (i,j) closes a multiloop.

In this case, the contribution to QB_i,j is given by the following (26) (27) (28)

Now QB_i,j = A_i,j + B_i,j + C_i,j. It nevertheless remains to define the recursions for QM1_i,j and QM_i,j. These satisfy the following.

(29)

(30)

(31)

This completes the derivation of the recursions for expected energy.

Comparison of structural entropy and derivational entropy

Using random RNA, 960 seed alignment sequences from Rfam family RF00005, and a collection of 2450 sequences obtained by selecting the first RNA from the seed alignment of each family from the Rfam 11.0 database [40], we show the following.

1. The thermodynamic structural entropy algorithms DP, FTD compute the same structural entropy values with the same efficiency, although as sequence length increases, FTD runs somewhat faster and returns slightly smaller values than does DP, since FTD uses a finite difference to approximate the derivative of the logarithm of the partition function.
2. DP and FTD appear to be an order of magnitude faster than the SCFG method of [24], which latter requires two minutes for RNA sequences of length 500 that require only a few seconds for DP and FTD.
3. Derivational entropy values computed by the method of [24] are much larger than thermodynamic structural entropy values of DP and FTD, ranging from about 4–8 times larger,depending on the SCFG chosen.
4. The length-normalized correlation between thermodynamic structural entropy values and derivational entropy values is poor to moderately weak.

Unless otherwise specified, throughout this paper, FTD, DP and SCFG refer to the formal temperature derivative method (Algorithm 1, with Turner’04 parameters), the dynamic programming method (Algorithm 2, with Turner’04 parameters), and the stochastic context free grammar method of [24]. SCFG(G4), SCFG(G5), SCFG(G6) respectively refer to the SCFG method of [24] using the stochastic context free grammars G4, G5, and G6. Additionally, there are three different training sets for each grammar: Rfam5, Mixed80 and Benchmark—see [24] for explanations of the training sets. Thus SCFG(G6,Benchmark) refers to derivational entropy, computed by the algorithm of [24], using grammar G6 with training set Benchmark, etc.

Table 1 lists the average values, plus or minus one standard deviation, for the entropy values and run time (in seconds) for 960 transfer RNAs from the seed alignment of family RF00005 from Rfam 11.0 [40]. Results for five methods are presented: (1) the dynamic programming method of this paper, using the Turner 2004 free energy parameters (DP), (2) approximating the formal temperature derivative by finite differences, and subsequently applying Eqs (10 and 6), using Turner 2004 free energy parameters (FTD); (3,4,5) using the program of [24] respectively with the stochastic context free grammars G4, G5, and G6 trained on the dataset ‘Rfam5’. For the 960 transfer RNAs from the Rfam database, this table shows that entropy values computed by DP and FTD are four to eight times smaller than derivational entropy values returned by the program of [24], while DP and FTD run five to ten times faster than the program of [24]—run times and derivational entropy values heavily depend on the grammar chosen and the training set used for production rule probabilities.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. Average values for structural entropy and run time (in seconds) for the 960 transfer RNA sequences from the seed alignment of Rfam family RF00005.

Methods include: DP: dynamic programming algorithm from our program RNAentropy, using the Turner 2004 energy parameters; FTD (ΔT = 10⁻⁷): finite difference computation of , where formal and table temperature are uncoupled, and formal temperature increment is 10⁻⁷; SCFG(G4,Rfam5): SCFG method [24] using grammar G4 with training dataset ‘Rfam5’ SCFG(G5,Rfam5): SCFG method using grammar G5 with training dataset ‘Rfam5’ SCFG(G6,Rfam5): SCFG method using grammar G6 with training dataset ‘Rfam5’. FTD returns very similar values for temperature increments 10⁻⁷ ≤ ΔT ≤ 10⁻¹¹; however, for smaller temperature increments, there is a slight deviation due to numerical precision issues—for example, average entropy of FTD with ΔT = 10⁻¹² is 5.238878±1.504748, with similar run times as other FTD runs.

https://doi.org/10.1371/journal.pone.0137859.t001

Table 2 presents the Pearson correlation for entropy values of 960 transfer RNAs from the seed alignment of family RF00005 from the database Rfam 11.0 [40]. The upper-triangular [resp. lower-triangular] entries are correlations for unnormalized [resp. length-normalized] entropy values. Entropy values were computed for the same methods as in Table 1. Since there is little variation in sequence length for the transfer RNAs in the seed alignment of RF00005 (average length is 73.41±5.13), any correlation due to sequence length is eliminated. The table shows the poor correlation between SCFG structural entropy, as computed by each grammar, with thermodynamic structural entropy.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 2. Pearson correlation for entropy values of 960 transfer RNAs from the seed alignment of family RF00005 from the database Rfam 11.0 [40].

Upper-triangular entries are for unnormalized entropy values, while lower-triangular entries are for length-normalized entropy values. Entropy values were computed for the same methods described in Fig 1; in particular, all SCFGs were trained with RF00005, as described in [24].

https://doi.org/10.1371/journal.pone.0137859.t002

Table 3 presents the average positional entropy, length-normalized structural entropy, and corresponding Z-scores for a small collection of experimentally confirmed conformational switches, collected by Giegerich et al. [41], and available on the RNAentropy web server. There appears to be no clear entropic signal for conformational switches, at least with respect to this small collection of sequences.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 3. Thermodynamic structural entropy, positional entropy, and corresponding Z-scores for a small collection of experimentally confirmed conformational switches, collected in [41]—sequences available at the RNAentropy web site.

For each sequence, the positional (resp. structural) entropy x was computed, along with the mean μ and standard deviation σ of 1000 dinucleotide shuffles of the sequence. The Z-score is then . Dinucleotide shuffles were computed, using the Altschul-Erikson algorithm [42] as implemented in [43]. Pearson correlation between Z-scores for positional and structural entropy is 0.4103.

https://doi.org/10.1371/journal.pone.0137859.t003

Table 4 presents the number of sequences, average length-normalized thermodynamic entropy, average entropy Z-score, average length-normalized ensemble defect, and average Z-score for for sequences in the seed alignment of several RNA familes from the Rfam 11.0 [40], as well as the precursor microRNAs from the repository MIRBASE [44]. Average values are given, plus or minus one standard deviation. The Z-score is defined as , where x is the entropy (resp. ensemble defect) of a given sequence, and μ (resp. σ denotes the mean (resp. standard deviation) of corresponding values for 100 random sequences having the same dinucleotides, obtained by using the Altschul-Erikson dinucleotide shuffling algorithm [42]. As shown by this table, Rfam family members appear to have lower structural entropy as well as ensemble defect than random RNA having the same dinucleotides, although the family RF00005 of transfer RNAs shows an exception to this rule for structural entropy. The most pronounced Z-scores for structural entropy and ensemble defect are for precursor microRNAs, which have very stable stem-loop structures. These results are generally comparable, with the exception of entropy Z-scores for RF00005, with results concerning minimum free energy (MFE) Z-scores from [43, 45]. Indeed, the particularly low MFE Z-scores of precursor miRNAs is used as a feature in the support vector machine miPred to detect microRNAs [46].

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 4. For several large families from the Rfam 11.0 database [40], and for MIRBASE precursor microRNA [44], the table presents the number of sequences (seq), length-normalized values of thermodynamic structural entropy (H) and ensemble defect (ens def), and the corresponding Z-scores for entropy and ensemble defect.

For each sequence from a given RNA family, 100 random sequences were generated with the same dinucleotides, using the Altschul-Erikson dinucleotide shuffling algorithm [42] as implemented in [43]—in the case of MIRBASE, only 10 random sequences were generated for each sequence. Subsequently, Z-scores were computed as , where x is the entropy (resp. ensemble defect) of a given sequence, and μ (resp. σ) is the mean (standard deviation) of 100 random sequences having the same dinucleotides.

https://doi.org/10.1371/journal.pone.0137859.t004

We now turn to the figures that support each of the four assertions made at the beginning of Section “Comparison of structural entropy and derivational entropy”. Fig 1 shows the average run times and entropy values for for DP, FTD (ΔT = 10⁻⁷), and the SCFG method of [24] using each of the grammars G4, G5 and G6 with training data from the set ‘Benchmark’. According to benchmarking work of [24] and [47], the grammar G6 seems somewhat better than G4 and G5. It is for this reason that we focus principally on the grammar G6, which was first introduced in the SCFG algorithm Pfold for RNA secondary structure prediction—see [27]. Fig 1A depicts average run times for DP, FTD, and SCFG methods, for 100 random RNA sequences of length n, where n ranges from 20 to 500 with an increment of 5. This figure shows that FTD and DP run faster by an order of magnitude than the SCFG methods—indeed, for length 500 RNAs, derivational entropy is computed in two minutes, while thermodynamic structural entropy is computed in a few seconds. The Fig 1B depicts the entropy values computed by DP, FTD (ΔT = 10⁻⁷), and SCFG methods. Note that for large RNA sequence length, entropy values returned by FTD are slightly smaller than those returned by DP, in agreement with the discussion in Section “Entropy by statistical physics”. Entropy values for the grammar G5 are considerably larger than those of FTD and DP, while entropy values for G4 and G6 are almost identical and approximately twice the size of those from G5.

Fig 2A presents graphs of length-normalized entropy values, computed by DP and SCFG. Using methods from algebraic combinatorics [48, 49], it is possible to prove that the length-normalized asymptotic structural entropy is constant, as observed in this figure. By numerical fitting, we find that the slope of the DP line is 0.087, while that of G6 is 0.329; i.e. SCFG entropy values using the G6 grammar are 3.78 times those of DP entropy. This is supported by Table 1, which suggests that G6 entropy values are 3.56 times larger than DP, while G4 and G5 entropy values are 6.71 resp. 6.85 times larger than DP entropy values. Fig 2B depicts the relative frequency of structural entropy values for DP, FTD, and SCFG methods for 960 transfer RNA sequences from the seed alignment of the Rfam 11.0 database [40].

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2.

(A) The average of length-normalized entropy values, as computed by DP and SCFG (G6,Benchmark), using the same data as described in the caption of Fig 1. Using methods from algebraic combinatorics, it can be proven that the length-normalized entropy for a homopolymer is asymptotically constant. By numerical fitting, we find that SCFG values are roughly four times as large as DP values (approximate fitted value 3.78). (B) Relative frequency of entropy values for the 960 transfer RNA sequences in the seed alignment of RF00005 family from Rfam 11.0 [40], as computed for each of the five methods DP, FTD (ΔT = 10⁻⁷), SCFG(G4,Rfam5), SCFG(G5,Rfam5) and SCFG(G6,Rfam5). See the caption from Fig 1 for explanation of each method, where in contrast to previous figures, the training set ‘Rfam5’ was used in place of ‘Benchmark’. Average entropy values for RF00005 are given as follows. FTD (ΔT = 10⁻⁷): 5.53±1.34. DP: 5.95±1.38; G4: 39.92±2.88; G5: 40.68±3.05; G6: 21.21±2.41. Note the bimodal distribution of entropy values computed with the SCFGs G4 and G5. Relative frequency plot for 712 5S ribosomal RNAs from RF00001 is very similar (data not shown).

https://doi.org/10.1371/journal.pone.0137859.g002

Fig 3 presents scatter plots and Pearson correlation of length-normalized entropy values and several notions of structural diversity that have been used for RNA design [14, 31]. Values were computed in this figure for a set of 2450 RNAs of various lengths, by selecting the first sequence from the seed alignment of each family from the Rfam 11.0 database [40], after discarding a few families having too few sequences. Fig 3A depicts the Pearson correlation between length-normalized structural entropy values, as computed by DP, FTD, and the SCFG method using grammars G4, G5, G6. Length-normalized derivational entropy values remain highly correlated, regardless of training set, but the correlation of all SCFG methods is poor with DP. The Pearson correlation of 0.79 for length-normalized entropy values obtained by G4 and G5 is high; however the correlation with G6 drops to 0.56 (G4-G6) and 0.34 (G5-G6). Fig 3B depicts scatter plots and Pearson correlation for 960 transfer RNAs from family RF00005 of Rfam 11.0, for length-normalized structural entropy, as computed by DP, and various notions of structural diversity used in synthetic RNA design. (By minimizing values such as the positional entropy, structural entropy, ensemble defect, expected base pair distance, it is more likely that computationally designed RNAs will fold into their predicted structures when experimentally validated.) Brief definitions of the notions of structural diversity that are compared in Fig 3B are given as follows. Native Contacts: proportion of base pairs in the Rfam consensus structure that appear in the low energy Boltzmann ensemble, defined by , where s₀ is the Rfam consensus structure. Positional Entropy: average positional entropy , where H₂(i) is defined by Eq (5). Expected base pair distance: length-normalized value determined from ∑_s p(s) ⋅ d _BP(s, s₀), where s₀ is the Rfam consensus structure, computed by ∑_{1 ≤ i < j ≤ n} I[(i, j) ∉ s₀] ⋅ p_{i, j} + I[(i, j) ∈ s₀] ⋅ (1−p_i,j) where I denotes the indicator function—see [14]. Ensemble defect: length-normalized value determined from , where s₀ is the Rfam consensus structure, I denotes the indicator function, and is defined in Eq (4). Vienna structural diversity: Boltzmann average base pair distance between each pair of structures in the ensemble, called ensemble diversity in the output of RNAfold -p[10], formally defined by ∑_{i < j} p_i,j(1−p_i,j)+(1−p_i,j)p_i,j, where p_i,j and output as ensemble diversity by RNAfold -p. Morgan-Higgs structural diversity: Boltzmann average Hamming distance between each pair of structures in the ensemble, where a structure s is represented by an array where s[i] = j if (i,j) or (j,i) is a base pair, and otherwise s[i] = i, formally defined by . Length-normalized DP entropy values are moderately highly correlated with positional entropy, but not with the other measures. In synthetic design of RNAs, it is our opinion that one should prioritize for experimental validation those synthetically designed RNAs by consideration of ensemble defect, structural entropy, etc., where the measures selected are not highly correlated. From this standpoint, one might use ensemble defect, structural entropy and proportion of native contacts as suitable measures for synthetic RNA design—see [50].

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3.

(A) Correlation between length-normalized structural entropy values, as computed by DP and five stochastic context free grammars: grammars G4, G5 and G6 for the ‘Benchmark’ training set, and G6 for ‘Rfam5’ and ‘Mixed80’ training sets (see [24]). Low correlation is shown between length-normalized thermodynamic structural and derivational entropies. For the fixed grammar G6, very high correlation is displayed between length-normalized entropy values for each of the training sets ‘Benchmark’, ‘Rfam5’, ‘Mixed80’ (similar results for fixed grammars G4,G5—data not shown). Although grammars G4 and G5 display a moderately high correlation together, there is low correlation with length-normalized entropy values determined by the grammar G6. Benchmarking set consists of the first sequence in the seed alignment from each family in the database Rfam 11.0 [40]. (B) Scatter plots and correlation between thermodynamic structural entropy and several measures of structural diversity, computed from 960 tRNA sequences in the seed alignment of family RF00005 from from the Rfam 11.0 database [40]. Correlation is computed between the following normalized values: (1) DP: length-normalized thermodynamic structural entropy computed by DP algorithm. (2) Native Contacts: proportion of base pairs in the Rfam consensus structure that appear in the low energy Boltzmann ensemble, defined by , where s₀ is the Rfam consensus structure. (3) Positional Entropy: average positional entropy, defined by , where H₂(i) is defined by Eq (5). (4) Expected base pair distance: length-normalized value determined from ∑_s p(s)⋅d_BP(s,s₀), where s₀ is the Rfam consensus structure, which equals ∑_{1 ≤ i < j ≤ n} I[(i,j) ∉ s₀]⋅p_i,j+I[(i,j) ∈ s₀]⋅(1−p_i,j) where I denotes the indicator function—see [14]. (5) Ensemble defect: length-normalized value determined from , where I denotes the indicator function, and is defined in Eq (4)—see [50]. (6) Str. Div. (V): Vienna structural diversity, output as ensemble diversity by RNAfold -p[10]. (7) Str. Div. (MH): Morgan-Higgs structural diversity [30], defined in the text. Positional entropy is moderately correlated with DP; ensemble defect and expected base pair distance are highly correlated, and each is moderately correlated with the proportion of native contacts. Structural diversity (Vienna and Morgan-Higgs) are highly correlated with positional entropy, but only (surprisingly) only moderately correlated with conformational entropy DP, in spite of the fact that all these measures concern properties of the ensemble of structures. Ensemble defect, expected base pair distance and expected number of native contacts are all highly correlated; this is unsurprising, since all measures concern the deviation of structures in the ensemble from the minimum free energy structure. Note that positional entropy is poorly correlated with the proportion of native contacts, although Huynen et al. [21] show that base pairs in the MFE structure of 16S rRNA tend to belong to the structure determined by comparative sequence analysis when the nucleotides have low positional entropy.

https://doi.org/10.1371/journal.pone.0137859.g003

Fig 4 displays the heat capacity and structural entropy for a thermoswitch (also called RNA thermometer) from the ROSE 3 family RF02523 from the Rfam 11.0 database [40], with EMBL accession code AEAZ 01000032.1/24229-24162. The heat capacity, computed by Vienna RNA Package RNAheat, presents two peaks, corresponding to two critical temperatures T₁,T₂, where one of the two conformations of this thermoswitch is stable in the temperature range between T₁ and T₂. The entropy plot also suggests the presence of a stable structure in the temperature range between T₁ and T₂, since small entropy values entail small diversity in the Boltzmann ensemble of structures.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. Heat capacity (left) and thermodynamic structural entropy (right) for a thermoswitch, or RNA thermometer, from the ROSE 3 family RF02523 from the Rfam 11.0 database [40], with EMBL accession code AEAZ01000032.1/24229-24162.

Lighter curves in the background correspond to the heat capacity (left) and thermodynamic structural entropy (right) of random RNAs having the same dinucleotides, obtained by the implementation in [43] of the Altschul-Erikson dinucleotide shuffle algorithm [42]. Since structural entropy and heat capacity the derivative of entropy H with respect to temperature closely follows the curve of the heat capacity (data not shown). Heat capacity computed using Vienna RNA Package RNAheat [10], and entropy computed by method DP.

https://doi.org/10.1371/journal.pone.0137859.g004

As shown in the tables and figures, the DP and FTD methods return almost identical values and have very similar (fast) run times, contrasted with the SCFG method, which is slow and whose values are much larger than those of DP and FTD. For a sequence of length 500, SCFG(G6,Benchmark) takes 2 minutes, compared with a few seconds for DP and FTD. Since FTD approximates a derivative by a finite difference, one expects a small discrepancy in the values of DP and FTD for thermodynamic structural entropy. According to [24], the sensitivity and specificity of G4 and G6 grammars are “significantly” higher than that of the G5 grammar. Since G6 is the underlying grammar of the Pfold software, for many of our comparisons, we compute derivational entropy using grammar G6 with the ‘Benchmark’ training set. (In data not shown, we benchmarked all nine combinations of grammars and training sets.)

Using RNAfold to compute conformational entropy

We have recently learned that newer versions of Vienna RNA Package [10] allow the user to modify the value RT by using the flag ––betaScale (kindly pointed out by Ivo Hofacker). It follows that RNAfold can easily be used to compute conformational entropy by using the FTD method. Let T = 310.15 be the absolute temperature corresponding to 37°C, let ΔT = 0.01, let T₂ = T + ΔT = 310.16 and T₁ = T − ΔT = 310.14. Define the scaling factors , and . Run RNAsubopt -p ––betaScale β₂ to compute the ensemble free energy −R(T + ΔT) ln Z(T + ΔT), and RNAsubopt -p ––betaScale β₁ to compute the ensemble free energy −R(T − ΔT) ln Z(T − ΔT), where Z(T + ΔT) [resp. Z(T − ΔT)] temporarily denotes the value of the partition function where table temperature is 37°C (as usual), and formal temperature is T + ΔT [resp. T − ΔT] in Kelvin. It follows that the uncentered finite difference Eq (32) (32) as well as the centered finite difference (33) both provide good approximations for the expected energy 〈E〉. Now run RNAsubopt -p to compute the ensemble free energy G = −RT ln Z where table and formal temperature are (as usual) 310.15 in Kelvin, and so compute the entropy (34)

Let ViennaRNA [resp. ViennaRNA*] denote the entropy computation just described, where expected energy is approximated by the uncentered Eq (32) [resp. centered Eq (33)]. Similarly, we let FTD [resp. FTD*] denote the uncentered [resp. centered] version of our code from Algorithm 1 in Section “Statistical Mechanics” in Methods. In computing entropy for Rfam family RF00005, both ViennaRNA and ViennaRNA* sometimes return entropy values that are larger than the correct values computed by DP, while entropy values of FTD [resp. FTD*] are always smaller than [essentially always smaller] than those of DP, as expected when using finite differences to approximate the derivative of the strictly decreasing, concave-down function ln Z(T).

Fig 5A shows that ViennaRNA is somewhat faster than FTD, and for each method, the uncentered version is faster than the centered version, which is clear since the former [resp. latter] computes the partition function twice [resp. three times]. Fig 5B shows that the standard deviation of entropy values for 100 random RNA is larger for ViennaRNA than FTD, and the uncentered form of ViennaRNA displays the largest standard deviation when ΔT = 10⁻⁴ (for ΔT = 0.01, all four finite derivative methods are comparable). These results are unsurprising due to numerical precision issues; e.g. for the 98 nt purine riboswitch with EMBL accession code AE005176.1/1159509-1159606, the algorithm DP determines a value of conformational entropy 9.975439, whereas by using (centered) ViennaRNA* with version 2.1.8 of RNAfold with ΔT = 10⁻², we obtain 9.93425742505. For ΔT = 10⁻⁴, 10⁻⁵ and 10⁻⁶, ViennaRNA* computes entropies of 9.59831636855, 6.94285165005⋅10⁻⁸, −5.9169597422⋅10⁻⁷. Such numerical instability issues are of much less concern to our method FTD and FTD*, as Fig 1 demonstrates for the uncentered method FTD with ΔT = 10⁻⁷.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 5. Average values for the run time and the entropy values for 100 random RNA sequences of length n, each having expected compositional frequency of 0:25 for A,C,G,U, where n ranges from 20 to 500 with increments of 5 for conformational entropy.

(A) Average run times as a function of sequence length, where error bars represent ±1 standard deviation. Methods used: DP, FTD, FTD*, ViennaRNA, ViennaRNA*. For random RNAs of length 500 nt, Vienna RNA Package is about three times faster than our code. (B) Standard deviation of the entropy values computed for 100 random RNA, displayed as a function of sequence length. From top to bottom, the first three curves represent uncentered ViennaRNA with ΔT = 10⁻⁴, centered ViennaRNA* with ΔT = 10⁻⁴, and DP. The bottom curve represents centered FTD with ΔT = 10⁻⁴, centered FTD* with ΔT = 10⁻², uncentered ViennaRNA with ΔT = 10⁻², centered ViennaRNA* with ΔT = 10⁻². The average entropy values computed by FTD, FTD*, ViennaRNA, and ViennaRNA* are indistinguishable and since FTD values are shown in the right panel of Fig 1, they are not shown here.

https://doi.org/10.1371/journal.pone.0137859.g005

Fig 6 shows the distribution of entropy differences (DP-FTD, DP-FTD*, DP-ViennaRNA, DP-ViennaRNA*) for 960 transfer RNAs from family RF00005 from the Rfam 11.0 database [40]. Reasons for the behavior of ViennaRNA and ViennaRNA* are presumably due to numerical precision issues. These differences are small, so when plotted as a function of sequence length in a manner analogous to Fig 1 (not shown), average entropy values computed by FTD, FTD*, ViennaRNA, and ViennaRNA* for ΔT = 10⁻² and 10⁻⁴ are visually indistinguishable.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 6.

(A) Relative frequency of the difference in entropy values for 960 transfer RNAs from the RF00005 family of the Rfam 11.0 database. (1) DP-FTD with average entropy difference 0.2512 ± 0.4935 with maximum of 3.1622 and minimum of 0. (2) DP-FTD* with average entropy difference 0.2502 ± 0.4934 with maximum of 3.1602 and minimum of -0.0020. (3) DP-ViennaRNA with average entropy difference 0.2475 ± 0.4975 with maximum of 3.1520 and minimum of -0.1743. (4) DP-ViennaRNA* with average entropy difference 0.2494 ± 0.4946 with maximum of 3.1572 and minimum of -0.0777. It is noteworthy that FTD is always less than DP, FTD* exceeds DP by a tiny margin only rarely, while ViennaRNA and ViennaRNA* more often exceed DP. Recall that the average deviation DP-FTD increases with increasing sequence length, as shown in the right panel of Fig 1. The same is true for DP-FTD*, DP-ViennaRNA, DP-ViennaRNA* (data not shown). (B) Free energy of arginyl-transfer RNA from Aeropyrum pernix with tRNAdb accession code tdbR00000589 [37] for temperatures ranging from 37°C to 38°C in increments of 0.01. The blue piecewise linear curve was created using RNAeval -T from the Vienna RNA Package [10]. The red linear curve was created by (1) calculating the entropy S_t = G(37) − G(38) of the tRNA cloverleaf structure by subtracting the free energy at 38°C from the free energy at 37°C, as determined using RNAeval -T, (2) computing the enthalpy H_t = G(37) + (273.15 + 37) · S_t, and then (3) computing the free energy at temperature T by G(T) = H_t − T · S_t. The jagged free energy curve is due to the fact that Vienna RNA Package represents energies as integers (multiples of 0.01 kcal/mol), so that loop energies jump at particular temperatures.

https://doi.org/10.1371/journal.pone.0137859.g006

Due to numerical stability issues, ViennaRNA and ViennaRNA* perform optimally with ΔT = 10⁻². Note that when using RNAfold, it is essential to use ––betaScale; indeed, if one attempts to compute the entropy using Eq (34) where expected energy is computed from Eq (32) [resp. Eq (33)] by running RNAfold -p -T 37.01 and RNAfold -p -T 37 [resp. RNAfold -p -T 37.01 and RNAfold -p -T 36.99], then the resulting entropy for the 98 nt purine riboswitch with EMBL accession code AE005176.1/1159509-1159606 is the impossible, negative value of -208.13 [resp. -210.61]. The large negative entropy values in this case are not only due to the lack of distinction between formal and table temperature, but as well to the fact that Vienna RNA Package represents energies as integers (multiples of 0.01 kcal/mol), so that loop energies jump at particular temperatures, as shown in the right panel of Fig 6. These issues should not be construed as shortcomings of the Vienna RNA Package, designed for great speed and high performance, but rather as a use of the program outside its intended parameters. As shown by Fig 5, the methods ViennaRNA and ViennaRNA* can rapidly compute accurate approximations of the conformational entropy.

Correlation with hammerhead cleavage activity

In [32], Shao et al. considered a 2-state thermodynamic model to describe the hybridization of hammerhead ribozymes to messenger RNA with subsequent cleavage at the mRNA GUC-cleavage site. In that paper, they define the total free energy (35) where each of these energies is defined on p. 10 of [32], and obtained by averaging over 1000 low energy structures sampled by Sfold[51]. The authors show a (negative) high correlation between ΔG_total and the cleavage activity of 13 hammerhead enzymes for GUC cleavage sites in ABCG2 messenger RNA (GenBank NM_004827.2) of H. sapiens; i.e. the lower the total change in free energy, the more active is the ribozyme. (Shao et al. originally considered 15 hammerheads; however two outlier hammerheads were removed from consideration.) Here, we show that the correlation with cleavage activity can be improved slightly by taking secondary structure conformational entropy into consideration.

To fix ideas, we consider the first GUC cleavage site considered by Shao et al. The minimum free energy (MFE) hybridization complex, as predicted by RNAcofold from the Vienna RNA Package [10] is shown in Fig 7A. The MFE structure of the 21 nt portion of mRNA, followed by a linker region of five adenines, followed by the hammerhead ribozyme, as computed by RNAfold from the Vienna RNA Package yields the same structure (where the linker region appears in a hairpin). It follows that to a first approximation, MFE hybridization structures can be predicted from MFE structure predictions of a chimeric sequence that includes a linker region. (Before the introduction of hybridization MFE software [10, 52], this approach was used to predict hybridization structures.)

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 7.

(A) Hybridization structure predicted by RNAcofold[59] of a 21 nt portion of messenger RNA for H. sapiens ABC transporter ABCG2 messenger RNA (GenBank. NM_004827.2) hybridized with a hammerhead ribozyme (data from the first line of Table 1 of [32]). The 21 nt portion of mRNA is 5′-UGCUUGGUGG UCUUGUUAAG U-3′ and the 42 nt hammerhead rizozyme is 5′-ACUUAACAAC UGAUGAGUCC GUGAGGACGA AACCACCAAG CA-3′. Messenger RNA is shown in green, while the hammerhead appears in red. In data not shown, we determined the secondary structure of the 21 nt mRNA portion, followed by a linker region of 5 adenines, followed by the 42 nt hammerhead ribozyme, by using RNAfold[10]. The base pairs in the hybridization complex are identical to the base pairs in the chimeric single-stranded sequence (not shown)—i.e. except for the unpaired adenines from the added linker region, the structures are identical. This fact permits us to approximate the structural entropy for the hybridization of two RNAs by using RNAentropy to compute the entropy of the concatenation of the sequences, separated by a linker region. (B) Correlation between hammerhead cleavage activity, as assayed by Shao et al. [32], with ΔG_d (change in free energy due to disruption of mRNA, denoted ΔG _disrupt in text), ΔG (change in total free energy, denoted ΔG _total in text), both taken from [32], with ΔS (change in conformational entropy k_B⋅ΔH), and ΔG(total) − T ΔS. Cleavage activity was measured by Shao et al. for the cleavage of GUC sites in ABC transporter ABCG2 messenger RNA (GenBank NM_004827.2). Values of ΔG_d, ΔG were taken from Table 1 of [32], while the change in conformational entropy ΔS was computed by RNAentropy. Note modest increase in the correlation of cleavage activity with ΔG, when adding the free energy contribution −TΔS, due to conformational entropy.

https://doi.org/10.1371/journal.pone.0137859.g007

In this case, enzyme activity is 0.843, ΔG _total = −5.423 kcal/mol, structural entropy of the hammerhead is 2.830, structural entropy of the 21 nt portion of mRNA is 2.146, and structural entropy of the 21 nt portion of mRNA portion with linker and hammerhead is 2.328. Assuming that the entropy of a rigid structure is zero, the change in structural entropy ΔH(hammerhead) is 0 − 2.830 = −2.830, and similarly ΔH(21 nt mRNA + linker) is −2.146, ΔH(21 nt mRNA + linker + hammerhead) is −2.328. The net change in structural entropy ΔH is ΔH(21 nt mRNA + linker + hammerhead) minus ΔH(21 nt mRNA + linker) minus ΔH(hammerhead), so ΔH = −2.328 − (−2.146 − 2.830) = 2.648. The net change in conformational entropy ΔS = k_B ⋅ ΔH is then 0.00526, hence the free energy contribution −TΔS = −RTΔH = −1.632. The correlation between ΔG _total and −TΔS is the value of 0.108, while the correlation value of −0.788 between hammerhead activity and ΔG _total is increase in absolute value to −0.806 (p-value 0.000878) when also taking into account −TΔS. See Fig 7 for a scatter plot and correlations between enzyme activity and ΔG [resp. ΔG − TΔS], which correspond to the total free energy change without [resp. with] a contribution from conformational entropy.

Fig 7A depicts the minimum free energy hybridization structure of a 21 nt portion of the ABC transporter ABCG2 messenger RNA from H. sapiens (GenBank NM_004827.2), hybridized with a hammerhead ribozyme (data from the first line of Table 1 of [32]). The MFE hybridization structure was computed by Vienna RNA Package RNAcofold[10]. We obtain the same structure by applying RNAfold to the chimeric sequence obtained by concatenating the 21 nt portion of mRNA, given by 5′-UGCUUGGUGG UCUUGUUAAG U-3′, with a 5 nt linker region consisting of adenines, with the 42 nt hammerhead rizozyme, given by 5′-ACUUAACAAC UGAUGAGUCC GUGAGGACGA AACCACCAAG CA-3′ (data not shown). By such concatenations with a separating 5 nt linker region, we can compute the structural entropy of hybridizations of the 21 nt mRNA with the hammerhead ribozyme. (In future work, we may extend RNAentropy to compute the entropy of hybridization complexes without using such linker regions.)

Structural entropy of HIV-1 genomic regions

Fig 8A depicts the structural entropy, computed as a moving average of 100 nt portions of the HIV-1 complete genome (GenBank AF033819.3). Using RNAentropy, the structural entropy was computed for each 100 nt portion of the HIV-1 genome, by increments of 10 nt; i.e. entropy was computed at genomic positions 1, 11, 21, etc. for 100 nt windows. To smooth the data, moving averages were computed over five successive windows. The figure displays the moving average entropy values, as a function of genome position (top dotted curve), entropy Z-scores, defined by , where x is the (moving window average) entropy at a genomic position, and μ [resp. σ] is the mean [resp. standard deviation] of the entropy for all computed 100 nt windows. Fig 8B is a portion of the NCBI graphics format presentation of GenBank file AF033819.3. Regions of low Z-score are position 4060 (Z-score of -2.69), position 8700 (Z-score of -2.46) and position 4040 (Z-score of -1.95). Since positions do not appear to correspond to the start/stop position of annotated genes, we ran cmscan from Infernal 1.1 software [53] on the HIV-1 genome (GenBank AF033819.3). We obtained 11 predicted noncoding elements as listed in Table 5, including the trans-activation response (TAR) element. Many of the predicted noncoding RNAs are much shorter than the 100 nt window used in the RNAentropy genome-scanning approach just described—it follows that low entropy Z-scores cannot be expected for such elements. Nevertheless, certain elements have quite low entropy Z-scores, such as the 5′-UTR and TAR element, both of which are known to be involved in the packaging of two copies of the HIV-1 genome in the viral capsid [54].

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 8. Structural entropy plot for the HIV-1 genome (GenBank AF033819.3).

UsingRNAentropy, the structural entropy was computed for each 100 nt portion of the HIV-1 genome, by increments of 10 nt; i.e. for 100 nt windows starting at genome position 1, 11, 21, etc. To smooth the curve, moving averages were computed over five successive windows. (A) Dotted-line displays moving average values of structural entropy; solid curve displays entropy Z-scores, defined by , where x represents the (moving window average) entropy at a genomic position, and μ [resp. σ] represents the mean [resp. standard deviation] of the entropy for 100 nt windows. Some of the lowest entropy Z-scores are -2.69 at position 4060, -2.46 at position 8700, -1.95 at position 4040. (B) NCBI graphics display of the HIV-1 genome, for comparison purposes. Low entropy (negative Z-score) regions do not appear to correspond with the start/stop location for annotated genes. In data not shown, we also computed positional entropy values [21] for the same windows, and determined a Pearson correlation of 0.7025 [resp. Spearman correlation of 0.6829] between (moving window average) values of entropy and positional entropy.

https://doi.org/10.1371/journal.pone.0137859.g008

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 5. Computationally annotated RNA noncoding elements from the HIV-1 genome with corresponding entropy Z-scores.

Running cmscan from Infernal 1.1 [53] on the HIV-1 genome (GenBank AF033819.3), we obtain 11 noncoding elements as listed in the table, along with the nucleotide beginning and ending positions, length of noncoding element, E-score, and entropy Z-score. Entropy Z-scores were computed using RNAentropy as explained in the text. Many of the annotated noncoding elements are much shorter than 100 nt, the length of the window size used; however, sporadic checking of entropy Z-scores computed for a moving window of size 50 does not seem to radically change the entropy Z-scores. Nevertheless, certain elements have low entropies and corresponding entropy Z-scores, such as the 5′-UTR and TAR (trans-activation response) element, both of which are known to be involved in the packaging of the HIV-1 genome in the viral capsid [54].

https://doi.org/10.1371/journal.pone.0137859.t005