A model of k-mer surprisal to quantify local sequence information content surrounding splice regions

Information theory

In 1948, Shannon linked the information content of a sequence of symbols, first described by Hartley, and entropy, a quantity used in thermodynamics (Shannon, 1948; Hartley, 1928; Gibbs, 1902). Shannon’s discoveries initially focussed on transmission of messages over noisy channels. These later became fundamental principles in information storage (MacKay, 2003). Shannon Entropy is a measure of the complexity of an ensemble X, of symbols x, where each symbol occurs with a probability p(x). The self-information associated with each symbol is called surprisal (Tribus, 1961) and is defined by: (1)${\displaystyle S\left(x\right)=-{log}_2\left(p\left(x\right)\right)}$ where S(x) is measured in bits. For the full ensemble of n symbols, the total information of the ensemble is the sum over all surprisals $I\left(X\right)={\sum }_{i=1}^{n}S\left(x_i\right)$. The Shannon entropy is defined as the average information per symbol or the expectation value of all surprisals: (2)${\displaystyle H\left(X\right)=E\left(S\left(x\right)\right)=-\sum _{i=1}^{n}p\left(x_i\right){log}_2\left(p\left(x_i\right)\right)}$ where ${\sum }_{i=1}^{n}p\left(x_i\right)=1$ and H(X) is also measured in bits. A few years after this formulation, the structure of DNA and the first protein sequences were discovered (Sanger, 1952; Watson & Crick, 1953). These discoveries and further advances in biology enabled the application of information theory to biological sequences. The storage of biological information within ensembles of molecular sequences led to several investigations of the quantification of biological information and its association with with biological functions (Gatlin, 1966). The Central Dogma of molecular biology itself was first framed in terms of the “transfer of sequential information” (Crick, 1970; Cobb, 2017), and information theory continues to underpin our understanding of how the genome encodes genetic information (Pritišanac et al., 2019).

A DNA sequence can be represented using an alphabet corresponding to the individual nucleotides 𝔸 = {A, C, G, T}. Since each amino acid is defined by a tri-nucleotide, or codon, it is useful to consider nucleotide sequences using an alphabet of 64 symbols, one for every possible tri-nucleotide. Of these, 61 codons encode amino acids, and the remaining three correspond to stop codons. Similarly, protein sequences can be defined by an alphabet of 20 symbols, one for each amino acid. This decline in the number of possible symbols between coding DNA, RNA, and protein leads to a decline in the maximum amount of information that can be encoded at each level. Information theory provides a theoretical framework within which to quantify these differences (Yockey, 1974; Nemzer, 2017).

The reason why almost all organisms translate only 20 amino acids despite the ability to encode 61 possible codons has not been fully determined (Koonin & Novozhilov, 2009). However, the redundancy in the genetic code allows for additional information to be captured within a coding region beyond that required to define the amino acid sequence itself (Yockey, 2000; Itzkovitz, Hodis & Segal, 2010). It has been suggested that codon degeneracy and the structure of the codon table support functions in addition to the encoding of amino acids (Dayhoff, Schwartz & Orcutt, 1978; Henikoff & Henikoff, 1992; Itzkovitz & Alon, 2007; Berleant et al., 2009; Maraia & Iben, 2014). These include the description of splicing regulatory motifs (Lim & Burge, 2001; Zhang & Chasin, 2004; Wang & Burge, 2008), DNA binding sites that co-exist within the coding sequence (Melnik & Usatenko, 2014; Vinga, 2014; Shreif, Striegel & Periwal, 2015), and RNA secondary structure (Itzkovitz, Hodis & Segal, 2010). This additional information can be viewed as a separate signal containing non-coding information multiplexed through the protein coding sequence.

While the genetic code naturally leads to a focus on triplet sequences, other representations are possible, and different length sequences reveal different aspects of the genome. For example, the information associated with individual nucleotides can be used to identify the presence of motifs within an ensemble of short molecular sequences (Schneider et al., 1986). It is important to consider that information content shown in motif figures is usually represented as 2 − H(X) as the aim is to identify consistent nucleotides within the motif rather than diversity (Schneider & Stephens, 1990). Dinucleotides can also be used to represent DNA sequences as they are important in the specification of epigenetic modifications (CpG islands), binding sites, splice donors and acceptors. Alphabets representing molecular sub-sequences of length k (k-mers) have also been widely used in motif discovery (Castle et al., 2008; Bailey, 2011). Information theoretic approaches have previously been used to investigate biological features such as coding and non-coding regions, nucleosome positioning and DNA binding sites with a variety of methods and representations to quantify information (Koslicki, 2011; Vinga, 2014; Wu, Zhang & Mu, 2014).

Splicing

Splicing is an essential mechanism in human cells performed by the spliceosome, a large ribonucleoprotein complex comprised of five small ribonucleoproteins at its core plus many other protein cofactors (Matera & Wang, 2014). The process of splicing involves precise removal of intron sequences from the transcribed RNA sequence. Each gene can express multiple mRNAs with different patterns of intron removal, such that exons in one mRNA may be part of an intron for another. Approximately 95% of genes are spliced in humans through this exquisitely regulated process, which is responsible for much of the diversity in the proteome (Lee & Rio, 2015; Sahebi et al., 2016). For a mechanistic review see Shi (2017). Exon splicing is determined by the binding of three key molecules: U1, splicing factor 1 (SF1) and the U2 auxiliary factor (U2AF) to the 5′ splice site, the branch point and 3′ splice site within the intron, respectively. The decision to include an exon within a transcript is generally made at the time of binding of these three molecules and is mediated by the use of Splicing Regulatory Elements (SREs)  (Sickmier et al., 2006; Diederichs et al., 2016; Saha et al., 2020). Some of these SREs are embedded within exons and can be viewed as loci where additional information is multiplexed along with the information required for amino acid sequence determination. Analysis of recurrent k-mers near splice junctions have already been shown to predict novel splice sites and sequences involved in splicing  (Lim & Burge, 2001; Zhang & Chasin, 2004; Fairbrother et al., 2004; Schwartz, Hall & Ast, 2009; Raponi et al., 2011; Ke et al., 2011; Erkelenz et al., 2014).

Here we have created a method whereby different biological sequences can be interrogated without requiring a direct multiple sequence alignment, or the need for sufficient sequence conservation with which to build that alignment. Using splice sites as a set of known biological features, we show that our model can quantify the information content of sequences at these regions irrespective of species.

Annotation

Genomic annotation for 5 different species Homo sapiens, Mus musculus, Drosophila melanogaster, Danio rerio, Schizosaccharomyces pombe was downloaded from Ensembl v99 (ftp://ftp.ensembl.org/pub/) (Cunningham et al., 2019). Genomic DNA, coding mRNA, and peptide files, were mapped against genomic annotation provided by the gene transfer format (.gtf) file. Only protein coding transcripts with GENCODE basic annotation were included in this analysis, which is defined as the set of 5′ and 3′ complete transcripts. Exon-intron boundaries and exon-exon junctions were independently generated using the gene transfer format annotation and sequences were aligned at splice sites.

Only k-mers with nucleotide symbols consisting fully of 𝔸_nt = {A, C, G, T} or 𝔸_aa = {A, C, D, E, F, G, H, I, K, L, M, N, P, Q, R, S, T, V, W, Y} for amino acids were used in this analysis. Since the majority of eukaryote genes have multiple transcripts per gene, with many exons shared between them, when evaluating coding sequences and amino acids, the k-mers from those shared exons would be duplicated if all transcripts were considered individually. We therefore impose a multiple transcript correction such that the genomic locations; chromosome, strand, start and end positions, for all k-mers must be unique. The repeated k-mers used by different transcripts at the same genomic locus are discarded. This method is also sensitive to k-mers spanning exon-exon junctions, since it maps the positions to genomic locations. The coding mRNA distribution contains all k-mers from within coding sequences only; introns, UTRs, and non-coding transcripts were excluded. Amino acid k-mers were reduced in length by a factor of 3, ( $k_{aa}=\frac{k_{nt}}{3}$), and their value was allocated to all three nucleotides in their corresponding codon when comparing with DNA and coding mRNA sequences.

Probability of k-mer occurrences

The k-mers starting at every position in the DNA (both strands), all coding mRNAs, and all amino acid sequences were extracted, counted and recorded in frequency tables. These frequency tables were then used to identify the probability of the k-mer, x, occurring in the overall sequence (3)${\displaystyle p\left(x\right)=\frac{counts\left(x\right)}{\sum _{x}counts\left(x\right)}}$ where ∑_xcounts(x) is equivalent to the total number of k-mers in the sequence, which corresponds to the sequence length minus k and any discarded sequences. Using this measure of k-mer probability, the surprisal for each k-mer can be calculated by Eq. (1) while the Shannon entropy for the total sequence can be found using Eq. (2).

Choice of k-mer length

This estimate of k-mer probabilities (Eq.(3)) works well for small k, however for larger k this entropy calculation is limited due to the finite sample size of biological sequences (Herzel, Ebeling & Schmitt, 1994; Herzel & Große, 1995). Since the total number of possible k-mers in nucleotide space is 4^k_nt and is 20^k_aa in amino acid space, the H. sapiens DNA sequence is too small to contain every possible k-mer for k_nt ≥ 17 in DNA sequences and k_aa ≥ 6 for amino-acids. The k-mer distributions associated with these spaces are also skewed, such that some k-mers do not occur even for much smaller k. Here we refer to k-mers with zero occurrences as ‘nullomers’ (Hampikian & Andersen, 2007). The presence of these nullomers at short k, combined with the skew in the distribution, lowers the observed entropy H^obs of the sequence. This information loss can be quantified in terms of the sequence redundancy: (4)${\displaystyle R=1-\frac{H^{obs}}{H^{max}}}$ where the largest possible entropy of the system, $H_{nt}^{max}=-{log}_2\left(4^{-k_{nt}}\right)=2k_{nt}$ bits for nucleotide sequences and $H_{aa}^{max}=-{log}_2\left(20^{-k_{aa}}\right)$ bits for amino acid sequences. Since H^max is dependent only on k whilst H^obs is limited by the length of the sequence, R increases as the total number of possible k-mers becomes greater than sequence size. This effect is also dependent on the uniformity of the k-mer distribution since H is maximised when all k-mers are equally likely.

Figure 1 shows the redundancy and the proportion of unique k-mers observed for 2 ≤ k ≤ 15 for the 5 different species: H. sapiens, M. musculus, D. rerio, D. melanogaster and S. pombe. As expected, all species show an increase in redundancy at larger k and a corresponding decrease in the number of k-mers represented at least once in the sequence. In DNA space (Figs. 1A and 1D), this effect becomes prominent at k_nt > 12. We therefore conclude that 12-mers are an appropriate size to model DNA sequences for all species except S. pombe, where they are too large relative to the size of the significantly smaller genome. We have therefore removed S. pombe from further analysis. Following similar reasoning, for comparisons between DNA, coding mRNA and amino acid sequences, we selected k_nt = 9 and k_aa = 3 to reflect the reduced amount of coding mRNA sequence relative to the DNA.

Surprisal patterns across splice sites

Metazoan exon-intron junctions have known conserved features within introns. Downstream of the 3′ end of the vast majority of exons is the GT donor, and upstream of the 5′ end of the exon is a corresponding AG acceptor site. A pyrimidine rich region occurs approximately 4–20 nt upstream of the acceptor site and a conserved but location variable branch point occurs approximately 15–55 nt upstream of the acceptor site (Corvelo et al., 2010). These motifs make exon-intron junctions ideal candidates for testing our model. Figure 2 shows the mean of DNA 12-mer surprisal across all protein coding exon-intron junctions for each species. A consistent surprisal pattern is observed, with substantial changes in the mean surprisal at the position of the polypyrimidine tract and splice site motifs. This result is expected since the surprisal, and hence information, is inversely proportional to k-mer occurrence, and therefore surprisal decreases in the presence of common sequences motifs. Conversely, Fig. 2 also reveals that for all species, exons contain significantly greater information than introns near splice boundaries. Although apparently intuitive, it is surprising since exon sequences are constrained to contain only those sequence patterns capable of representing a functional protein, while intron sequences are under no such constraint. For example, k-mers featuring an in-frame stop codon are not permissible within a coding sequence. A naïve view of coding space, therefore, is one in which exon sequences are constructed from a subset of possible k-mers while intron sequences can be constructed from the entire repertoire of k-mers. This predicts that exons would be constructed from more common sequences, and that their surprisal would therefore be lower than that of the corresponding intronic space. By contrast, these data suggest a greater sequence constraint on intronic regions near splice junctions than similar proximal exon regions. This in turn implies that there are additional constraints for intronic sequences, such as those arising from the need to encode intronic SREs. It is tempting to speculate that these patterns result from a selection pressure that excludes certain valid coding sequences from intronic space. Importantly these patterns are robust against different values of k, demonstrating the generality of the model (Fig. S1). Here we show data for k = 12, since these have a lower coefficient of variation than those of smaller k (Fig. S2).

Surprisal across exon junctions

Since the distinct patterns in Fig. 2 arise from differences between coding and non-coding sequence spaces, we aimed to identify additional information encoded in the DNA and coding sequences beyond that which is required simply to define the encoded protein sequence. Protein coding exon-exon junctions were aligned at 9-mer resolution (amino acid 3-mers), with k = 9 chosen to account for the smaller sample size of coding sequences (Fig. 3). The mean surprisal for H. sapiens surrounding exon-exon junctions is 18.0 bits in DNA and 17.4 bits for coding mRNA, which as expected, is slightly lower due to the loss of in-frame stop codon sequences. Amino acid sequences show a significant reduction of approximately 1/3 information when compared to DNA sequences (12.6 bits). This is in keeping with previous work that considers the entropy of codons (Yockey, 1974). All other species considered here are consistent with H. sapiens, and show similar patterns in DNA surprisal across exon-exon junctions, with a constant difference dependent on the species. All species show a sharp decrease in surprisal at position −10 and sharp increase at position −1, 1nt upstream of the exon junction, consistent with Fig. 2. These 9-mers are created by the juxtaposition of the last nucleotide of the 5′ exon with the first 8 nucleotides of the 3′ exon. An interesting observation is that the most frequent dinucleotide across the junction for both exon-exon and exon-intron boundaries is GG, suggesting that the joining of the G in position −1 with the following 8-mer is driving this peak. For coding mRNA most species show a tight and consistent pattern similar to that of DNA surprisal, however D. melanogaster is a clear outlier, with a different pattern and a significantly higher information content. This is likely to be in part a consequence of the increased proportion of coding sequence in D. melanogaster (22.1% versus < 3% for the other species in Fig. 3).

The sharp decrease in surprisal is at position −9 for coding mRNA surprisal. When computed at other values of k (Fig. S3), these minima shift with k, indicating that the nucleotides driving this pattern are at the end of the k-mer, and suggesting that this pattern arises from the lesser conserved exonic 3 ′ AG terminal motif. There is little variation surrounding splice junctions for amino acid sequences which is expected since the information content encoded at these positions is expected to be a feature of splicing information, which is ahead of the translation process. This is consistent for other values of k (Fig. S3) and similarly to Fig. S2, the coefficient of variation is reduced for larger k (Fig. S4). Together these data strongly suggest that codon redundancy allows additional signals to exist within the DNA without significant impact on the amino acid sequence encoded through the same space.

Effect of CG dinucleotides

The spectra of k-mers in the DNA sequences of several species were previously described by Chor et al. (2009). Figure 4 shows the H. sapiens 12-mer DNA spectrum, which reveals three peaks in k-mer abundance corresponding to 12-mers occurring 4, 27, and 298 times in the genome. The spectrum can be viewed as three overlapping distributions, which are largely explained by the number of CG dinucleotides that occur within each 12-mer. This behaviour is not observed for any other dinucleotides (Fig. S5).

This CG dependency was also reported by Chor et al. (2009) for k_nt = 8, who suggested that it is a property of tetrapod genomes as a consequence of a global repression of CG dinucleotides, as shown in Fig. 4C for H. sapiens. This effect is further exacerbated by the concentration of CG dinucleotides at CpG islands, represented by a small number of CG-rich, common k-mers. The repression of CG dinucleotides is an established feature of genome evolution and results from the high mutation rate of CG dinucleotides caused by the deamination of methylated cytosines (reviewed by Walsh & Xu (2006)). A secondary consequence of the disproportionate rarity of CG dinucleotides is that k-mers that contain them tend to have high surprisal (Fig. 4D).

The frequency of CG dinucleotides in proximity to exon-intron junctions was investigated for H. sapiens DNA sequences (Fig. 5). Similar to Fig. 2, the ‘valley’ observed in mean surprisal centred at −16 is due to the effect of the polypyrimidine tract, shown explicitly in Fig. 5C. The CG dinucleotide frequency is 2-fold enriched within exons, providing a partial explanation for the increased information content in exons compared with introns. These results are in keeping with previous reports describing high CG differential between exon and intron regions (Amit et al., 2012). However, even within exonic regions, the CG dinucleotide occurs at a much lower frequency than all other dinucleotides. No other dinucleotides show such a disproportionate difference between intron and exonic spaces (Fig. S6). However, CG dinucleotide usage is not sufficient to explain all patterns, as shown by the lack of correspondence between CG dinucleotide usage and surprisal across the polypyrimidine tract.