Introduction

Studies of genomes based on linguistic approaches date a few decades back (Brendel et al. 1986; Pevzner et al. 1989; Searls 1992; Botstein and Cherry 1997; Gimona 2006; Faltýnek et al. 2019; Ji 2020). An interplay with methods of statistical physics as well as theory of complex systems brought new insights into biology (Dehmer et al. 2009; Qian 2013). Studies range from attempted n-gram-based classification of genomes (Tomović et al. 2006; Huang and Yu 2016) to algorithms for optimal segmentation of RNAs in secondary structure predictions (Licon et al. 2010) and analysis of substitution rates of coding genes during evolution (Lin et al. 2019), just to mention a few. Various types of sequences in genomes are related to multiple genetic codes (Trifonov et al. 2012) and can be studied both using quantitative linguistic point of view (Ferrer-i-Cancho et al. 2013; Ferrer-i-Cancho et al. 2014) and from a wider perspective, within more abstract approaches (Neuman and Nave 2008; Barbieri 2012). Recently, neural networks and deep learning algorithms emerged as new tools to analyze nucleotide sequences (Fang et al. 2019; Singh et al. 2019; Melkus et al. 2020; Ren et al. 2020) offering wider prospects for studies of genomes. Viruses, balancing on the fuzzy border between non-alive and alive, hence remaining on the verge of life (Villarreal 2004; Kolb 2007; Carsetti 2020), are within the most interesting subjects of studies.

The aim of the present Letter is to draw attention to simple treatments of nucleotide sequences in viral RNAs by means of new parameters, which can be immediately extracted from genome data. We expect that such parameters can be potentially used as an auxiliary tool in the classification of viruses (cf., in particular, Wang 2013). The idea of this study is linked to the recent COVID-19 outbreak, and the analysis started from comparing human coronaviruses (Su et al. 2016; Wu et al. 2020) and some other viruses. To achieve relative homogeneity of the material, we restrict our sample to single-stranded RNA viruses only. Both positive- and negative-sense RNAs are considered. For future reference, we also include two retroviruses, HIV-1 and HIV-2.

The paper is organized as follows. Summary of data and description of methods are given in “Data and Methods” section. Results are presented in “Results” section. Finally, brief discussion is given in “Discussion” section.

Data and Methods

The viral genomes are taken from the databases of the National Center for Biotechnology Information (NCBI, https://www.ncbi.nlm.nih.gov); the complete list is given in Table 1. Note that coronaviruses have rather long RNA genomes of ca. 30 kilobases (kba), which might bias the values of calculated parameters. To study the effect of RNA sizes, we also include some very short genomes, namely, Hepatitis D virus with 1682 ba (Saldanha et al. 1990) and Phage MS2 virus with 3569 ba (de Smit and van Duin 1993), as well as two longest known RNA viruses, Ball python nidovirus with 33452 ba (Gorbalenya et al. 2006) and Planidovirus with 41178 ba (Saberi et al. 2018). Still, the sizes of RNA viruses are much more homogeneous (the difference is up to 25 times) than those of DNA ones, which may vary by about four orders of magnitude (Campillo-Balderas et al. 2015).

Table 1 Viruses analyzed in the work
Full size table

We use two approaches to define nucleotide sequences. The first one is based on cutting an RNA genome into chunks of equal length of n nucleotides. The second approach is rooted in linguistics, so that the most frequent nucleotide is treated as a “space” dividing a RNA into “words” of different lengths (Rovenchak 2018). Note also distantly related units applied in the analysis of the human DNA, so called motifs (Liang 2014).

To demonstrate the first approach, with equal-length chunks, let us consider the Ebolavirus genome, starting with the following nucleotide sequence:

$$ \text{GGACACACAAAAAGAAAGAAGAATTTTTAGGATCTTTTGT\ldots}\ . $$
(1)

Choosing the chunk length n = 4, we obtain:

$$ \mathrm{GGAC\ ACAC\ AAAA\ AGAA\ AGAA\ GAAT\ TTTT\ AGGA\ TCTT\ TTGT\ \ldots}\ . $$
(2)

Eventually, for RNA length not being multiples of four, the last chunk can have one to three nucleotides. Obviously, the number of all possible 4-nucleotide combinations is 44 = 256. Note that longer chunks would yield much higher variety of combinations with frequencies being distributed very smoothly. On the other hand, we would like to avoid studies of shorter chunks, like three-nucleotide sequences corresponding to codons. So, the length n = 4 seems optimal for our analysis.

In the second approach, the same Ebolavirus sequence (1) can be split using the most frequent nucleotide – adenine – as a “space” into the following:

$$ \mathrm{GG\ C\ C\ C\ X\ X\ X\ X\ G\ X\ X\ G\ X\ G\ X\ TTTTT\ GG\ TCTTTTGT\ \ldots}\ . $$
(3)

The “X” stands for a zero-length element inserted between two consecutive “A”s.

We have also applied peculiar treatment of the Influenza A virus (H1N1) by adding spaces between each of eight segments of its RNA in the first and second approaches.

In both approaches, we calculate the frequencies of obtained nucleotide chunks within a given genome split in the respective manner and compile the rank–frequency distributions. The latter are obtained in a standard manner as follows: the most frequent item has rank 1, the second most frequent one has rank 2 and so on. Items with equal frequencies are given consecutive ranks in a random order, which is not relevant.

Results

The rank–frequency distributions obtained using the first approach – with 4-nucleotide chunks – were analyzed using a special software, AltmannFitter 2.1 (Altmann 2000). We found that two discrete distributions describe the obtained data with the highest precision, so called 1-displaced negative hypergeometric distribution (Grzybek 2007; Wilson 2013):

$$ p_{r}=\frac{ \left( \begin{array}{cc}{M+r-2}\\{r-1}\end{array}\right) \left( \begin{array}{cc}{K-M+n-r-2}\\{n-r+1}\end{array}\right) }{ \left( \begin{array}{cc}{K+n-1}\\{n}\end{array}\right) }, \qquad r=1,2,3,\ldots\ . $$
(4)

and Pólya distribution (Wimmer and Altmann 1999; Johnson et al. 2005):

$$ p_{r}=\frac{\left( \begin{array}{cc}{-p/s}\\ {r-1}\end{array}\right) \left( \begin{array}{cc}{(p-1)/s}\\ {n-r+1}\end{array}\right) }{ \left( \begin{array}{cc}{-1/s}\\ {n}\end{array}\right) }, \qquad r=1,2,3,\ldots\ . $$
(5)

Absolute frequencies are obtained by multiplying pr by the sample size N. In most cases, the discrepancy coefficient C = χ2/N is smaller than 0.02, which is considered a good fit (Mačutek 2008). Typical rank–frequency distributions and respective fits are shown in Fig. 1. Complete data are summarized in Table 2 and visualized in Fig. 2.

Fig. 1
figure 1

Typical rank–frequency distributions of four-nucleotide chunks and respective fits. The left panel shows the data for MERS and the fit with the hypergometric distribution, which is one of the best (C = 0.0011). The right panel demonstrated the worst fit obtained for the Hepatitis D virus data fit with the Pólya distribution (C = 0.0342)

Full size image
Table 2 Fitting parameters for the distributions of four-nucleotide chunks
Full size table
Fig. 2
figure 2

Location of viruses on the K − M plane (negative hypergeometric fit, left panel) and s − p plane (Pólya fit, right panel)

Full size image

Note that the abovementioned goodness-of-fit condition comes from quantitative linguistic domain, where it is an “empirical rule of thumb” (Antić et al. 2019), and might not be directly applicable in the studies of genomes. However, it is used for various language levels, including letters and words (Antić et al. 2019; Mačutek 2008), and the observed distribution of four-nucleotide chunks are very similar to those of letters or syllables (cf, Wilson 2013; Rovenchak et al. 2018).

The first immediate observation from Fig. 2 is that the length of genomes has no special influence on the fitting parameters. Indeed, both the shortest Hepatitis D genome and two longest – Ball python nidovirus and Planidovirus – genomes have close values of M or s parameters. On the other hand, for genomes of similar lengths (coronaviruses) a clear separation is seen with respect to M and p parameters. It is even more pronounced in the former case corresponding to the negative hypergeometric distribution: lower values for HCoV viruses (229E, HKU1, NL63, and OC43) and higher ones for MERS, SARS, and SARS-CoV-2.

Rank–frequency distributions were also compiled for nucleotide “words” obtained using the second approach and used to calculate certain parameters, like entropy, mean length (first central moment), length dispersion (second central moment) and some others. A typical rank–frequency distribution is shown in Fig. 3. Comparing it to Fig. 1 we can easily see qualitatively different behavior, in particular, significantly longer plateaus at high ranks / low frequencies, which makes such distributions close to those of words in human languages. The differences between the distribution shapes of nucleotide “words” and n-grams are pronounced especially well in samples of rather short sizes, like several kba or a few dozen kba, and thus properties of “word-like” sequences might give new data for studies of such material, including mitochondrial DNA and viral RNA.

Fig. 3
figure 3

Typical rank–frequency distribution of nucleotide “words” (the data for the Marburg virus are shown). The left panel shows the plot in the log–log scale while in the right panel linear scales over axes are used

Full size image

Previous studies (Rovenchak 2018) showed that entropy and mean lengths of such “word-like” nucleotide sequences in the mitochondrial DNA can be used to distinguish species and genera of mammals. It appears, however, that even better results are achieved with the “entropy – length dispersion” pair of variables, cf. Fig. 4.

Fig. 4
figure 4

Grouping of mammal species on the m2 − S plane. Red-shaded area corresponds to Felidae, the blue one denotes Ursidae, and the green-one corresponds to Hominidae. Calculations are made using mitochondrial DNAs with adenine as a “space”

Full size image

The parameters are defined as follows. Entropy is given by

$$ S=-\sum\limits_{r=1}^{r_{\max}} p_{r} \ln p_{r} , $$
(6)

where the upper summation limit corresponds to the total number of different “words” in the list and relative frequencies pr are

$$ p_{r}=f_{r}/N,\qquad\text{where}\ N=\sum\limits_{r} f_{r} $$
(7)

and fr are absolute frequencies at rank r. Mean length and length dispersion are

$$ m_{1}=\frac{1}{N} \sum\limits_{i} x_{i}, \qquad m_{2}=\frac{1}{N} \sum\limits_{i} (x_{i}-m_{1})^{2}. $$
(8)

where the summations run over all the “words” of the analyzed genome. Lengths xi of a particular word are counted as the number of nucleotides except for “X” having length zero.

One should note that from similarity of species one can expect proximity of points but not vice versa: it would be too bold to expect species distinguishability from only two parameters.

This second approach can be divided into two sub-branches: (a) adenine, which is the most frequent nucleotide in most species studied in the present work, is used as a “space”; (b) the most frequent nucleotide is used as a “space”. The latter is mostly relevant for RNAs, where low frequencies of adenine yield too long “words” thus significantly distorting the expected dependencies. The respective results are shown in Figs. 5–7. All the data are summarized in Table 3.

Table 3 Parameters for the distributions of nucleotide sequences separated by a specific nucleotide
Full size table

In Fig. 7, we can observe in particular that α-coronaviruses, HCoV-229E and HCoV-NL63, have very close values of the parameters (the respective point nearly overlap). A similar situation is with β-corovaniruses HCoV-OC43 and HCoV-HKU1. Two other β-corovaniruses, SARS and SARS-CoV-2, are located close to HCoV-OC43 and HCoV-HKU1, while MERS occupies an intermediate position. The latter virus also significantly differs in the entropy value, see Fig. 5. On the other hand, calculations with the most frequent nucleotide used as a space (T for the analyzed coronaviruses) do not exhibit such a grouping, see Fig. 6.

Fig. 5
figure 5

Location of viruses on the m2 − S plane. Calculations are made using RNAs with adenine as a “space”, hence entropy is denoted SA

Full size image
Fig. 6
figure 6

Location of viruses on the m2 − S plane. Calculations are made using RNAs with the most frequent nucleotide as a “space”, hence entropy is denoted Sm.f.

Full size image
Fig. 7
figure 7

Location of viruses on the m2 − S/N plane. Calculations are made using RNAs with adenine as a “space”, hence entropy is denoted SA. The vertical axis thus represents the entropy divided by the number of nucleotide sequences separated by adenine in the respective genome

Full size image

Similarly to the case of fixed-length chunks (four-nucleotide sequences analyzed above), one can expect close points for similar species but should not deduce that close points mean related species. Figures 4–6 demonstrate only one pair of parameters obtainable from rank-frequency distributions of nucleotide “words”, while Table 3 contains additional data. Further analysis can be done by processing the complete raw dataset used for calculations, which is freely available at https://doi.org/10.5281/zenodo.4045875.

When looking in detail into the rank–frequency distributions corresponding to coronaviruses we have discovered the following pattern: the first rank is always occupied by “X” followed by three single-nucleotide “words” with ranks 2–4, while the fifth ranks are occupied by a two-nucleotide sequence with either the same (4-same) or different (4-diff) nucleotides, see Table 4. Curiously, different nucleotides correspond to coronaviruses causing much more severe diseases. This observation is yet to be extended onto a wider material, but the preliminary data for the analyzed human viruses are as follows:

  • 4-same: Dengue, HCoV-229E, HCoV-HKU1, HCoV-NL63, HCoV-OC43, HIV-1, HIV-2, HRV-A, HRV-B, HRV-C, Polio;

  • 4-diff: A/H1N1, Ebola, Hepatitis A, Hepatitis C, Hepatitis E, Marburg, Measles, MERS, Norovirus, Rabies, SARS, SARS-CoV-2.

Three other viruses, Hepatitis D, Yellow fever, and Zika, do not follow either pattern having a two-nucleotide sequence with as low ranks as 3 or 4.

Table 4 Top-ranked nucleotide sequences in the genomes of the human coronaviruses
Full size table

The abovementioned feature can be viewed with respect to (dis)similarity in the frequency structure of viral RNA and RNA or DNA of infected species. For instance, the fifth most frequent sequence of human mtDNA split by the most frequent C nucleotide is TA. Interestingly, it coincides with a sequence composed of nucleotides paring with G and C – and the GC sequence is the fifth most frequent in, e.g., MERS, SARS, SARS-CoV-2, see Table 4.

Discussion

We have presented several possible approaches to simple parametrization of RNA viruses based on the analysis of nucleotide sequences in viral genomes. They are based on discrete distributions (negative hypergeometric and Pólya) for equal-length (4-nucleotide) chunks and on the pair “entropy – length dispersion” for distributions of sequences separated by adenine or another most frequent nucleotide. Close values of parameters calculated from rank-frequency distributions of various nucleotide sequences are characteristic for related viruses, which is connected to similar structures of viruses and thus might reflect similarities in their functional properties. In some cases, however, close values are also obtained for unrelated viruses. This is not surprising as representing viruses on a plane means a two-parametric projection of points that are certainly described by more than two variables. We consider our study as preliminary steps in discovering such variables.

The structure of rank–frequency distributions of sequences of both types suggest an analogy with two text levels. Namely, equal-chunk sequences approach to letters or syllables while nucleotide “words” obtained using certain nucleotide as a separator are similar to ordinary words. Here, we should stress that these data correspond to rather short molecules, viral RNA and previously studied mtRNA (Rovenchak 2018), so the correlations might differ when dealing with larger nucleotide strings (cf. Gimona 2006; Ferrer-i-Cancho et al. 2013; Faltýnek 2019). Such an interpretation corresponds to the biosemiotic analogy between natural language texts and strings of biopolymers.

Observations regarding peculiarities of rank–frequency distributions, with the fifth most frequent sequence containing two either the same or different nucleotides (4-same vs 4-diff), support the fact that 4-diff cases correspond to viruses causing potentially more severe diseases when dealing with seven human coronaviruses. This tendency is generally preserved if the analyzed set is expanded by other viruses studied in this work. Some precautions concern, in particular, the two HIV types, which fall into the 4-same category while certainly being extremely dangerous. However, HIV are not strictly RNA viruses but retroviruses, so we suggest that the reported peculiarities might be specific for RNA viruses only. “False-positive” alerts (cf. Norovirus in the 4-diff category) are not problematic, but the rate of “false-negative” results (severe diseases in the 4-same category) is yet to be identified. Expansion of the analyzed material in future studies would help to clarify the relevance of this observation. To establish relations between peculiarities of the rank–frequency distributions in virus genomes and disease severity, a formalization of the latter is required. Initially we planned using the case fatality rate (CFR) indicator (Reich et al. 2012; Kim et al. 2020) but where not able to find a study with data for different viruses based on a unified approach, similar, e.g., to (GBD 2017).

The main expected outcome of our reported analysis is a call for collaboration to expand the dataset and consistently classify diseases caused by RNA viruses, in particular with respect to severity and contagiousness. If some simple patterns could be established in the nucleotide distributions, this might help alerting healthcare systems, which seems to become a very topical issue from this year on.