RESEARCH ARTICLE

Frequency Analysis Techniques for Identification of Viral Genetic

Data
Vladimir Trifonov and Raul Rabadan
Department of Biomedical Informatics, Center for Computational Biology and Bioinformatics, College of Physicians and Surgeons, Columbia University, New York, New
York, USA

ABSTRACT Environmental metagenomic samples and samples obtained as an attempt to identify a pathogen associated with the
emergence of a novel infectious disease are important sources of novel microorganisms. The low costs and high throughput of
sequencing technologies are expected to allow for the genetic material in those samples to be sequenced and the genomes of the

Downloaded from http://mbio.asm.org/ on December 7, 2018 by guest

novel microorganisms to be identified by alignment to those in a database of known genomes. Yet, for various biological and
technical reasons, such alignment might not always be possible. We investigate a frequency analysis technique which on one
hand allows for the identification of genetic material without relying on alignment and on the other hand makes possible the
discovery of nonoverlapping contigs from the same organism. The technique is based on obtaining signatures of the genetic data
and defining a distance/similarity measure between signatures. More precisely, the signatures of the genetic data are the frequen-
cies of k-mers occurring in them, with k being a natural number. We considered an entropy-based distance between signatures,
similar to the Kullback-Leibler distance in information theory, and investigated its ability to categorize negative-sense single-
stranded RNA (ssRNA) viral genetic data. Our conclusion is that in this viral context, the technique provides a viable way of dis-
covering genetic relationships without relying on alignment. We envision that our approach will be applicable to other microbial
genetic contexts, e.g., other types of viruses, and will be an important tool in the discovery of novel microorganisms.
IMPORTANCE Multiple factors contribute to the emergence of novel infectious diseases. Implementation of effective measures
against such diseases relies on the rapid identification of novel pathogens. Another important source of novel microorganisms is
environmental metagenomic samples. The low costs and high throughput of sequencing technologies provide a method for the
identification of novel microorganisms by sequence alignment. There are several obstacles to this method, as follows: our knowl-
edge of biology is biased by an anthropomorphic view, microbial genomic material could be a minuscule fraction of the sample,
the sequencing and enrichment technologies can be a source of errors and biases, and finally, microbes have high diversity and
high evolutionary rates. As a result, novel microorganisms could have very low genetic similarity to already known genomes, and
the identification by alignment could be computationally prohibitive. We investigate a frequency analysis technique which al-
lows for the identification of novel genetic material without relying on alignment.

Received 2 June 2010 Accepted 20 July 2010 Published 24 August 2010

Citation Trifonov, V., and R. Rabadan. 2010. Frequency analysis techniques for identification of viral genetic data. mBio 1(3):e00156-10. doi:10.1128/mBio.00156-10.
Editor Ian Lipkin, Columbia University
Copyright © 2010 Trifonov and Rabadan. This is an open-access article distributed under the terms of the Creative Commons Attribution-Noncommercial-Share Alike 3.0
Unported License, which permits unrestricted noncommercial use, distribution, and reproduction in any medium, provided the original author and source are credited.
Address correspondence to Vladimir Trifonov, [email protected].

N ext-generation high-throughput sequencing technologies

have significantly increased our ability to generate genetic
data at low costs. The expectations are that costs will decrease
Tool from NCBI) or SHRiMP (SHort Read Mapping Package from
Computational Biology Lab, University of Toronto) to compare the
genetic material with that in an existing reference database. The suc-
while the throughput and quality levels increase. The availability cess of this approach depends on the existence of a database contain-
of such data has allowed us to make significant progress in under- ing a reference with a sufficiently high level of similarity. For example,
standing the genetic bases of biological processes and systems (1). if the sample is obtained from a host species, this approach allows for
Metagenomic studies, such as those resulting from the use of en- identification of the host genetic material, if a reference for the host is
vironmental samples (2) or the human biome (3), are one particular available.
direction of research that was enabled by the new sequencing tech- This paper focuses on techniques towards categorization of
nologies. An important biochemical and bioinformatic challenge genetic material in the absence of a high-homology reference. In
presented by such studies is categorization of the genetic material this context, genetic material of viral origin presents a particular
present in the sample according to its species of origin. In the context challenge. At present, the NCBI database contains hundreds of
of identifying novel pathogens (4, 5), the species of origin can be a thousands of viral genomes. Despite the amount of effort put into
distant evolutionary relative to a known pathogen, and identification building such genomic databases, it is unlikely that they are even
of this relative is an important step towards understanding the patho- close to being comprehensive. In fact, further extension of these
gen present in the sample. A first step towards such categorization is databases is one of the goals of metagenomic studies.
using alignment tools like BLAST (Basic Local Alignment Search The situation is complicated further from three directions. On

July/August 2010 Volume 1 Issue 3 e00156-10 mbio.asm.org 1

Trifonov and Rabadan.

one hand, the reference database is biased towards known human have references with very low homology due to biased sampling
pathogens. One can argue that this bias is inherent to the way the and high mutation rates at our disposal. One cannot help but
samples are obtained— often as a result of an effort to fight the compare this task to deciphering ancient scripts, such as the
disease associated with the pathogen. Although metagenomic Mayan, Ugaritic, and Sumerian scripts, etc., without relying on a
samples are often obtained to study pathogen activity, the chal- close reference. The success of these archeological and linguistic
lenging cases occur when the pathogen cannot be identified and achievements can only be an inspiration for us. For the rest of this
isolated by traditional biochemical methods, which means that it paper, we present the techniques which were successfully applied
is unlikely that it will be included in the reference database. Fur- to the recent identification of a novel viral pathogen affecting
thermore, in some metagenomic studies, samples without any no- farmed salmon (7), although they were developed for dealing with
ticeable pathogen activity are obtained. In such cases, the chance this seemingly daunting task in a general setting.
of matching novel viral species present in the sample to known
references is diminished further. One can hope that as the cost of RESULTS
obtaining sequenced metagenomic data decreases, the influence Even though RNA viruses have relatively small genomes at about

Downloaded from http://mbio.asm.org/ on December 7, 2018 by guest

of such sampling biases will decrease as well. 104 nucleotides, the number of possible sequences is enormous,
The second obstacle goes deeper because it stands in the way of many more than the number of quarks in the universe. The di-
even extracting the viral genetic material from a sample. The issue mension of this space is significantly reduced due to correlations
is that the total amount of viral genetic material present in a sam- imposed by selective pressures, but even accounting for this di-
ple is dwarfed by other genetic material, host or otherwise, also mension reduction, the diversity of the viral world is high. An idea
present. The genetic material of an RNA virus consists of about of the level of diversity can be gained by using the generalized
104 to 105 nucleotides, which is 3 to 4 orders of magnitude less Shannon entropy, as defined in reference 8. This entropy in the
than a typical prokaryotic genome and 4 to 5 orders less than a case of the hemagglutinin gene of avian influenza A virus is 0.60,
typical eukaryotic genome. Current high-throughput sequencing and for the PB2 polymerase gene of human influenza A virus, this
technologies proceed roughly by first fracturing the genetic mate- entropy is 0.08, a result of 90 years of evolution. In contrast, in the
rial present in a sample, amplifying the fractured sample by PCR, context of human evolution, the entropy in the P53 gene is 0.21, a
and sequencing the PCR product. If the microbe is not grown in result of more than 100 million years of evolution.
culture and only clinical samples are available, the genetic material The problem with categorizing viral genetic material is related
from the microbe is usually at very low concentrations, and the to the problem of assigning a similarity or distance measure on
available sequences do not cover its whole genome. Biases and this space. The goals one has in mind when designing such a met-
errors introduced by the sequencing and enrichment processes are ric come from two sometimes contradictory directions: on the one
other problems which contribute to an incomplete picture of the hand, we want the metric to capture as much evolutionary con-
viral genetic material present in a sample. Part of such biases is the nectedness as possible; on the other, we want the measure to be
fact that the high throughput of the new sequencing technologies efficiently computable. One natural measure is the edit distance—
comes at the price of producing reads of very short lengths (6). the number of nucleotide mutations (substitutions, insertions,
This, on one hand, limits the genome reconstructing ability of and deletions) necessary to transform the genome of one virus
alignment techniques and, on the other, restricts the power of into another. This metric is at the core of traditional alignment
techniques, such as ours, based on the statistical properties of the algorithms, e.g., the Needleman-Wunsch and Smith-Waterman
reads. algorithms.
While the previous obstacles can be surmounted purely by The NCBI BLAST tool is one popular implementation of an
technological advancement, the third difficulty in the way of cat- alignment algorithm. Efficiency considerations in this algorithm
egorizing viral genetic data is far more fundamental because it is have led to the requirement that the genomes being compared
inherent to viruses themselves. Viruses depend on their hosts for contain a common seed of some predetermined length, which for
reproduction, and often, the survival of a viral pathogen is at the randomly chosen genomes corresponds to high similarity. For
expense of its host. The strongest weapon in the viral arsenal example, a seed of 10 nucleotides corresponds to 90% genetic
against the highly developed and sophisticated defense mecha- similarity in randomly chosen genomes. Considering that the al-
nism of the host is an extremely high mutation rate. It is estimated gorithm employed by BLAST represents more or less the top of the
that the mutation rate of a typical RNA virus is 10⫺4 nucleotides line in efficient alignment techniques for comparing genomes, one
per replication. If a virus takes a few hours to reproduce and it can say that low nucleotide homologies due to more distant evo-
makes thousands of offspring, in a year the descendant viruses can lutionary relatedness for such techniques are out of reach (9). To
differ substantially from their ancestors. One can gain perspective overcome this problem, a number of algorithms based on the
on the amount of genetic variability existing in the viral world by search of characteristic oligonucleotide motifs and profiles were
considering the evolution of the H1N1 influenza A virus from developed (10–12).
1918 to that of the recent 2009 H1N1 pandemic. The variability in The general framework of our approach is to project the high-
this case is comparable to the genetic distance between the human dimensional genetic space to a low-dimensional numeric space
and mouse genomes (~15%), a result of 100 million years of evo- and then define a similarity measure in this space. We refer to the
lution. The implication is that even if a related reference is present low-dimensional projection as a signature. For the purpose of this
in a viral database, the genetic distance to this reference can be paper, the signature of a given genome is its k-mer content, with k
prohibitively large for capturing the similarity by alignment. being a natural number. Thus, the signature of a genome is a
To summarize, in order to categorize viral genetic material 4k-dimensional vector. This approach to comparison was first ex-
present in sequenced metagenomic data, we have to identify the plored by Blaisdell (13), and the work of Vinga and Almeida (14)
species of origin for sequences with low coverage and for which we contains a review of algorithms based on it. For this study, we

2 mbio.asm.org July/August 2010 Volume 1 Issue 3 e00156-10

 Frequency Analysis for Identification of Viral Data

Downloaded from http://mbio.asm.org/ on December 7, 2018 by guest

FIG 1 The performance of FASD with various parameters. The top row contains the fraction of queries with the correct target in a certain percentile, the middle
row is the same but with the percentile of the top-scoring member of the same family, and the bottom row contains queries and targets that are disjoint.

concentrated on cases when k equals 1, 2, or 3. The main reason for mine the target sequences closest to the query for some fixed dis-
this is that viral genomes are very short, on the order of 10,000 tance measure on the signature space. The algorithm simply or-
nucleotides for RNA viruses. Furthermore, due to low concentra- ders the target sequences from the database according to their
tions of viral material and the resulting low coverage obtained distance to the query. The Jenson-Shannon divergence was ex-
from current sequencing technologies, candidate viral contigs are plored previously as a possible distance measure between genomic
very short, at lengths of approximately 500 to 1,000 nucleotides. In signatures (16). In this work, we concentrate on a symmetrized
principle, our technique is extensible to larger k if sufficiently long version of the Kullback-Leibler distance, defined below in Mate-
contigs are available. The guideline one should keep in mind is rials and Methods, which is motivated by likelihood consider-
that to have good estimates of the frequencies of k-mers from a ations. Unless specified otherwise, we refer to this notion of dis-
genome of length L, L must be at least several times larger than the tance as the frequency distance. The Kullback-Leibler divergence,
4k value. as a measure of similarity among genetic sequences, was also stud-
We investigated several possible distance measures to assess ied by Wu et al. (17).
their ability to categorize viral genetic data. First, maximum like- To analyze the performance of frequency analysis of sequence
lihood considerations led us to the entropy-based distances data (FASD), we took all the reference negative-sense single-
(Kullback-Leibler [KL], Jensen-Shannon, and ␭ divergences) of stranded RNA (ssRNA) viral segments available at NCBI (279 se-
the frequency estimates obtained from the signatures. Second, as quences) to form our target sequence database. For a fixed length,
Karlin and Burge (15) point out, since the relative dinucleotide from each sequence in the target database we extracted 100 equally
abundance is characteristic of the species producing the genetic spaced subsequences of such lengths to obtain a database of que-
material, we considered the rectilinear distance of the relative ries with 27,900 sequences. For every percentile, we obtained the
dinucleotide abundances, as defined by those authors (15). Third, fraction of query sequences that had their target sequences of or-
we considered the Euclidean distance and the correlation coeffi- igin ranked in that percentile in the order produced by FASD. We
cients of the frequency vectors because they form natural mea- performed this analysis for 1-, 2-, and 3-mers with queries having
sures of similarity with well-established general mathematical lengths of 500 nucleotides (Fig. 1, top left) and for 3-mers with
properties. Fourth, we considered the ␹2 test statistic as a measure queries having lengths of 300, 500, and 800 nucleotides (Fig. 1, top
of the likelihood of the distribution of k-mers in the two genomes right). This analysis was repeated with respect to the highest per-
happening by random chance only. centile of a viral segment from the same family as the sequence of
Given a query nucleotide sequence and a database of target origin of the query (Fig. 1, middle). In the case of 3-mers with
nucleotide sequences, the goal of frequency analysis is to deter- queries having lengths of 500 nucleotides, we found that for 88%

July/August 2010 Volume 1 Issue 3 e00156-10 mbio.asm.org 3

Trifonov and Rabadan.

of the queries, the correct target is in the top 5%, and for 97% of gets that come from different segments (Fig. 2, bottom left). Sec-
the queries, the correct target family is in the top 5%. ond, we compared the distribution of distances for all pairs of
In the previous analysis, every query sequence is part of a target queries coming from the same segment to the distribution of dis-
sequence. To analyze the performance of FASD without the as- tances of queries from different segments (Fig. 2, bottom right).
sumption that the database contains a target to which the query The sensitivity and specificity of a simple test, which takes a
aligns, we formed a disjoint target database in which for every threshold frequency distance and decides whether two sequences
query subsequence, we included the remaining, complementary are similar if they are closer than this threshold, can be assessed
part of the sequence of origin of the query. The disjoint target with the relative operating characteristic (ROC) curve of the test.
database again contains a correct target sequence for every query Based on the distribution of distances given in Fig. 2, we computed
sequence, but now the query and its correct target are disjoint and the ROC curves of this test for the following two cases: (i) when the
so do not align. We performed the same percentile analysis as positive examples are pairs of negative-sense ssRNA viral se-
described in the previous paragraph, except that for every query, quences and the negative examples are pairs of a negative-sense
we considered only the list of targets to which it did not align well. ssRNA viral sequence and an rRNA sequence and (ii) when the

Downloaded from http://mbio.asm.org/ on December 7, 2018 by guest

This list necessarily contained the correct target sequence, because positive examples are nonaligning pairs from the query and a
it was disjoint from the query. To determine whether two se- complementary target from the same negative-sense ssRNA virus
quences aligned well, we used the Smith-Waterman local align- and the negative examples are nonaligning pairs from the query
ment algorithm, with the following scoring: match, 2; mismatch, and complementary target from different negative-sense ssRNA
⫺3; gap opening, ⫺5; and gap extension, ⫺3. The choice of align- viruses (see the description of the disjoint target database above).
ment parameters was influenced by the existing NCBI computed The ROC curves in both cases are far from random tests, whose
parameters for obtaining the corresponding Karlin-Altschul ROC curves will follow the diagonal line. By optimizing specificity
E values. Our main goal was to catch almost exact alignments, and sensitivity, we find that in the first case, at a frequency distance
hence the high penalties for gap opening and extension. For a of 0.23, we achieve a false-positive rate of 22% and a true-positive
given alignment score, we obtained its Karlin-Altschul E value rate of 77% (Fisher’s exact test P value of ⬍10⫺7), and in the
(18), with parameters for the computation of the E value obtained second case, at a frequency distance of 0.19, we achieve a false-
from NCBI. We considered two sequences to align well if the positive rate of 32% and a true-positive rate of 69% (Fisher’s exact
E value was at most 10⫺2. The results of this analysis for 3-mers test P value of ⬍10⫺8).
with lengths of 300, 500, and 800 nucleotides are shown in Fig. 1, Another tactic we used to analyze the performance of FASD
bottom left. For 3-mers with queries at lengths of 500 nucleotides, was to mutate randomly the sequences of the target database. Our
we also compared the performances of FASD in the two target model of mutation has only one parameter—the expected frac-
databases (Fig. 1, bottom right). As expected, in the disjoint case, tion of mutated sites. For a given fraction m, mutation of a se-
the performance of FASD degrades, and in particular for 3-mers quence involves going over the sites of the sequence and indepen-
with queries having lengths of 500 nucleotides, we found that for dently deciding with probability m whether the site is mutated.
60% of the queries, the correct target is in the top 5%, as opposed For every site chosen for mutation, we pick a new nucleotide uni-
to 88% of those in the overlapping case. Nevertheless, even in the formly from the three available choices. For a given mutation pa-
disjoint case, the dependence of the cumulative fraction of queries rameter, we can form a mutated target database where we ran-
on the percentile of the correct target exhibited in the bottom row domly mutate each sequence in the target database according to
of Fig. 1 is far from trivial and stays close to the dependence when that parameter. For 3-mers with queries at lengths of 500 nucleo-
there is alignment (Fig. 1, bottom right). tides, we performed the percentile analysis in mutated target da-
We also investigated the distribution of the frequency distances tabases with 0%, 20%, 30%, and 40% mutated sites and computed
between negative-sense ssRNA viruses and from negative-sense the distribution of distances for 3-mers at lengths of 500 nucleo-
ssRNA viruses to other sets of genetic data. Given two sets of tides (Fig. 3, top).
sequences, we computed the pairwise frequency distance between Finally, we compared the performance of FASD with the fre-
pairs of sequences, with one from the first set and the other from quency distance to its performance with the ␹2 test statistic, the cor-
the second set. For one of the sets, we fixed the query database of relation coefficient, the Euclidean distance, and the Karlin-Burge dis-
27,900 of sequences with lengths of 500 nucleotides obtained from tance defined in reference 15 and the Materials and Methods section.
negative-sense ssRNA viruses as explained above. For the other In each case, we subtracted the corresponding cumulative distribu-
set, we first took the disjoint targets database and considered only tion of the percentile of the correct target from the same distribution
pairs including a query and a target, which do not align. Second, computed using the frequency distance (Fig. 3, bottom). As can be
we took the large-subunit (LSU) rRNA sequences included in the seen, the frequency distance slightly outperforms the other distance
Silva database and, finally, the segments of positive-sense ssRNA measures in this setting.
viruses deposited in NCBI (600 sequences). In each case, we com-
puted the distribution of 3-mer frequency distances (Fig. 2, top DISCUSSION
left). Frequency analytic techniques provide an approach to sequence
To understand how pieces of the same virus relate to each other similarity, which does not rely on alignment. Alignment tech-
(as opposed to pieces of the different viruses), we compared the niques, such as the one used in BLAST, have become standard
distribution of distances of pieces from the same virus to that of tools for sequence comparison when the similarity is high enough.
pieces from different viruses. More precisely, we first computed However, the comparison between two genomes is problematic
the distribution of distances between all pairs of queries and their when the relationship is distant. This is a problem that appears
corresponding complementary targets and compared it to the dis- often in the identification of novel/emergent pathogens using
tribution of distances between all queries and complementary tar- high-throughput sequencing technologies. High evolutionary

4 mbio.asm.org July/August 2010 Volume 1 Issue 3 e00156-10

 Frequency Analysis for Identification of Viral Data

Downloaded from http://mbio.asm.org/ on December 7, 2018 by guest

FIG 2 The top left row contains the distribution of frequency distances between negative-sense ssRNA viruses and other sources of genetic material. The top
right row contains the ROC curves for the comparison of distributions of distances for negative-sense ssRNA viruses versus rRNA and same versus different
segments; the best trade-off between FPR (false positive rate ⫽ number of false positives/total number of negatives) and TPR (true positive rate ⫽ number of true
positives/total number of positives) is marked. The bottom row shows the distribution of distances between pairs of subsequences from the same virus versus that
from pairs of subsequences from different viruses.

rates, lack of knowledge about related genomes, and the high rate An application of the method is to relate different genes or
of sequence errors of current sequencing technologies complicate segments from the same virus, even when they do not align at all.
the identification of the few reads/contigs from the pathogen ge- For instance, by applying FASD to the PB1 gene of the 2009 H1N1
nome that could be present in a clinical sample. pandemic, we can identify not only the close relatives of this gene
In this work, we studied a frequency analysis method for com- in several influenza A viruses but also other segments from the
parison of genomic data that is not based on sequence alignment same viruses as well (Table 1). It is equivalent to using BLAST to
but on the statistical properties of the genetic sequences repre- align hemagglutinin and also finding neuraminidase.
sented by their signatures. These properties take into account Another application of the method is to “assemble” reads/con-
codon biases, mutational biases, e.g., errors in polymerase activity, tigs from the same organism without necessarily requiring an
and selection biases, e.g., pressure from the immune system of the overlap between them. We refer to this form of alignment as “hor-
host (19, 20) and restriction sites in phages infecting a bacterial
izontal” to distinguish it from the “vertical” alignment required
host. In addition, the frequency analysis method is robust under
for the traditional assembly algorithms. As an example of horizon-
frame shift sequencing errors (e.g., the homopolymer errors in-
tal assembly, we took two viruses, parainfluenza and rabies, and
troduced by pyrosequencing) and genomic rearrangements.
split their genomes into 10 disconnected nonoverlapping pieces.
In the viral context, since viruses depend on their host for replica-
tion, pressures from the host are an important source of bias in viral The tree in Fig. 4 gives the result of a hierarchical clustering using
genetic signatures. These pressures are a result of a complex virus- the average method of combining clusters performed on the pair-
host interaction, which includes the particulars of the replicative and wise frequency distances of the pieces. It contains two distinct
infectious mechanism of the virus, e.g., the cellular context of the clusters corresponding to the two viruses of the example.
interaction, the structure of the virus, and its evolutionary history, Another example of the ability of FASD to relate sequences that
including previous and other coexisting hosts. A better understand- cannot be aligned is given by the plot in Fig. 5. Here we have
ing of the factors contributing to signature biases will provide a better computed the pairwise frequency distances and the logarithms of
understanding of the boundaries within which the frequency method the Karlin-Altshul E values of the Smith-Waterman alignment
succeeds. In this work, we assume the existence of such biases and scores between the coding sequences of the influenza A, B, and C
leverage them towards categorization of viral genetic data. viruses and the Ebola Zaire virus. As can be seen, the frequency

July/August 2010 Volume 1 Issue 3 e00156-10 mbio.asm.org 5

Trifonov and Rabadan.

Downloaded from http://mbio.asm.org/ on December 7, 2018 by guest

FIG 3 The top row shows the performance of FASD with respect to various levels of mutation (mut) according to a simple model of mutation. The bottom row
compares the performance of FASD with frequency distance and other distance measures of frequency vectors. corr, correlation coefficient.

distance can often indicate evolutionary relatedness and/or simi- analysis performed well and was able to detect genetic relation-
lar sources of origin, although the sequences do not align. ships even in the cases where there was no alignment. Due to its
Frequency analysis techniques, based on the k-mer composi- nature, the frequency analysis approach is less specific than an
tion of genomic data, provide an alternative way of uncovering alignment-based one. In fact, one can imagine the two approaches
biological relationships, different from standard techniques based lying on the opposite sides of a spectrum, with frequency analysis
on sequence alignment. The technique is applicable to cases in relying on common statistical properties of short subsequences
which the genetic data allows for good estimates of the k-mer and alignment relying on the presence of a common long subse-
frequencies, e.g., when contigs and genomes are sufficiently long. quence. The lack of specificity of frequency analysis is both its
We found that in the simulated setting of our study, the frequency weakness, as it produces false positives, and its strength, as it can

TABLE 1 Top 13 negative-sense ssRNA viral sequences produced by FASD, given as a query of the PB1 gene of the 2009 H1N1 pandemic
Frequency Value for L ⫻
distance frequency distance P valuea Viral segment
0.0011 5.12 3.14E⫺01 A/New York/392/2004 (H3N2) segment 2
0.0028 12.86 3.05E⫺02 A/Korea/426/68 (H2N2) segment 2
0.0041 16.41 9.86E⫺03 A/goose/Guangdong/1/96 (H5N1) segment 4
0.0047 21.73 1.76E⫺03 A/Hong Kong/1073/99 (H9N2) segment 2
0.0051 20.53 2.61E⫺03 Influenza B virus RNA 5
0.0056 25.75 4.73E⫺04 A/goose/Guangdong/1/96 (H5N1) segment 2
0.0059 27.21 2.92E⫺04 A/Puerto Rico/8/34 (H1N1) segment 2
0.0061 24.59 6.92E⫺04 A/Puerto Rico/8/34 (H1N1) segment 4
0.0062 28.44 1.95E⫺04 A/Hong Kong/1073/99 (H9N2) segment 1
0.0071 32.32 5.38E⫺05 A/Korea/426/68 (H2N2) segment 1
0.0073 29.19 1.52E⫺04 A/Korea/426/68 (H2N2) segment 4
0.0076 33.69 3.41E⫺05 A/New York/392/2004 (H3N2) segment 3
0.0077 28.27 2.06E⫺04 Maize fine streak virus, complete genome
a The computation of the P value is explained in “Statistics” in Materials and Methods.

6 mbio.asm.org July/August 2010 Volume 1 Issue 3 e00156-10

 Frequency Analysis for Identification of Viral Data

nucleotides indexed by all k-mers, such that CX is the number of

times the k-mer X appears in the sequence S.
Estimation of frequencies. Given a genomic signature C of a
sequence S, the vector F of frequencies of k-mers appearing in S is
obtained first by adding one to each of the components of C to
obtain a vector P of pseudo-counts. Then, letting L be the sum of
the components of P, the frequency of the k-mer X is calculated as
follows: FX ⫽ PX/L.
Frequency distance between genomic signatures. The frequency
divergence between two nucleotide sequences, Sa and Sb, is the
Kullback-Leibler divergence between the frequency vectors Fa and Fb

div(Sa, Sb) ⫽ 冘F
X a,X log2 Fa,X ⁄ Fb,X

Downloaded from http://mbio.asm.org/ on December 7, 2018 by guest

The frequency divergence is always nonnegative and is equal
to zero if and only if Fa equals Fb. Its usage is motivated by
FIG 4 Result of hierarchical clustering using the average method on the
considering the likelihood of sequence Sa occurring by chance
pairwise 3-mer frequency distances of 10 nonoverlapping genomic subse-
quences at lengths of 500 nucleotides from two viruses (1 to 5 from rabies virus if we pick a random sequence with k-mer frequencies fixed like
and 6 to 10 from parainfluenza 1 virus). those of sequence Sb. The frequency distance between se-
quences Sa and Sb is given by the following version of the
Kullback-Leibler distance:
uncover deeper genetic connections. In some cases, frequency
dist(Sa, Sb) ⫽ [La div(Sa, Sb) ⫹ Lb div(Sb, Sa)] ⁄ (La ⫹ Lb)
analysis could be the only viable way of detecting genetic relation-
ships in data obtained from high-throughput sequencing experi- where La and Lb are the sums of the components of pseudo-count
ments. We believe that frequency analysis is an important tool in vectors of Sa and Sb, correspondingly. The frequency distance is
the fight against emergent infectious diseases and the discovery of symmetrical, always nonnegative, and equal to zero if and only if
novel microorganisms. Fa equals Fb. The frequency distance does not satisfy the triangle
inequality.
MATERIALS AND METHODS Karlin-Burge distance. For a given sequence, let F1 and F2 be
Genomic signatures. For a given nucleotide sequence S and a its 1-mer and 2-mer frequency vectors. For nucleotides X and Y,
natural number k, the signature of S is a vector C of length 4k let the equation KXY ⫽ F2,XY/(F1,X·F1,Y) be the relative 2-mer fre-

FIG 5 The left plot contains the frequency distances between pairs of segments from four viruses (FLUA, influenza A; FLUB, influenza B; FLUC, influenza C;
EBOVZ, Ebola virus Zaire). The right plot is the logarithms of the Karlin-Altschul E values computed for the Smith-Waterman alignment scores between those
pairs of segments.

July/August 2010 Volume 1 Issue 3 e00156-10 mbio.asm.org 7

Trifonov and Rabadan.

quency of the 2-mer XY compared to its frequency as predicted by We are especially grateful to the researchers at the Center for Infection
the frequencies of the 1-mers X and Y. Then, the Karlin distance and Immunity for many interesting discussions and insights. We also
for two sequences, Sa and Sb, with relative 2-mer frequencies of Ka want to thank Miguel Brown for helping us to create a web interface.
and Kb, respectively, is calculated as follows:
distK(Sa, Sb) ⫽ (1 ⁄ 16) 冘 XY |Ka,XY ⫺ Kb,XY|
REFERENCES
1. Mardi, E. 2008. The impact of next-generation sequencing technology on
Statistics. Frequency analysis relies on the Kullback-Leibler genetics. Trends Genet. 24:133–141.
(KL) divergence as a measure of distance between a query and a 2. Williamson, S., D. B. Rusch, S. Yooseph, A. L. Halpern, K. B. Heidel-
target genome. By thinking of the KL divergence as a statistic of a berg, J. I. Glass, C. Andrews-Pfannkoch, D. Fadrosh, C. S. Miller, G.
Sutton, M. Frazier, and J. C. Venter. 2008. The sorcerer II global ocean
random event, we need to assess its significance. We will show that sampling expedition: metagenomic characterization of viruses within
for sufficiently long random query genomes, the distribution of aquatic microbial samples. PLoS One 3:e1456.
this statistic is approximated by a gamma distribution. 3. Gill, S., M. Pop, R. T. Deboy, P. B. Eckburg, P. J. Turnbaugh, B. S.
Take a query genome of length L with a k-mer frequency vector Samuel, J. I. Gordon, D. A. Relman, C. M. Fraser-Liggett, and K. E.

Downloaded from http://mbio.asm.org/ on December 7, 2018 by guest

q and a target genome with a k-mer frequency vector t. Assuming Nelson. 2006. Metagenomic analysis of the human distal gut microbiome.
Science 312:1355–1359.
that each k-mer of the query genome is chosen independently 4. Lipkin, I., G. Palacios, and T. Briese. 2007. Emerging tools for microbial
from the rest with probability as in t, the probability of obtaining diagnosis, surveillance, and discovery, p. 413– 435. In S. M. Lemon, et al.
a genome with k-mer counts like those of the query genome fol- (ed.), Global infectious disease surveillance and detection: assessing the
lows a multinomial distribution. Let M equal 4k. For large k-mer challenges—finding solutions. National Academies Press, Washington,
counts, using Sterling’s approximation to the factorial function, DC.
5. Travassos da Rosa, A. P., T. N. Mather, T. Takeda, C. A. Whitehouse,
this probability can be approximated by
R. E. Shope, V. L. Popov, H. Guzman, L. Coffey, T. P. Araujo, and R. B.
1⫺M Tesh. 2002. Two new rhabdoviruses (Rhabdoviridae) isolated from birds
(2␲L) 2 ⫺L · K
写
qi e
i
during surveillance for arboviral encephalitis, northeastern United States.
Emerg. Infect. Dis. 8:614 – 618.
6. Pop, M., and S. Salzberg. 2008. Bioinformatics challenges of new se-
where K is the KL divergence of the distribution q from the distri- quencing technology. Trends Genet. 24:142–149.
7. Palacios, G., M. Lovoll, T. Tengs, M. Hornig, S. Hutchison, J. Hui, R. T.
bution t. The probability of obtaining a genome of length L with Kongtorp, N. Savji, A. V. Bussetti, A. Solovyov, A. B. Kristoffersen, C.
KL divergence at most the value K is approximated by the sum of Celone, C. Street, V. Trifonov, D. L. Hirschberg, R. Rabadan, M.
the probabilities of all such genomes. For large lengths (L), this Egholm, E. Rimstad, and W. I. Lipkin. 2010. Heart and skeletal muscle
sum is close to the corresponding integral, and this integral can be inflammation of farmed salmon is associated with infection with a novel
approximated with the following observations. The KL ball of ra- reovirus. PLoS One 5:e11487.
8. Ricotta, C., and L. Szeidl. 2006. Towards a unifying approach to diversity
dius K can be approximated by an ellipsoid volume with radii measures: bridging the gap between the Shannon entropy and Rao’s qua-
(2K)1/2(t11/2. . . tM1/2) centered at t. Since we are integrating over dratic index. Theor. Popul. Biol. 70:237–243.
frequency vectors, this M-dimensional volume is intersected by a 9. Koski, L., and G. Golding. 2001. The closest BLAST hit is often not the
hyperplane through its center, resulting in an integral over an (M nearest neighbor. J. Mol. Evol. 52:540 –542
⫺ 1)-dimensional ellipsoid volume. Appropriate coordinate 10. Teeling, H., A. Meyerdierks, M. Bauer, R. Amann, and F. Glöckner.
2004. Application of tetranucleotide frequencies for the assignment of
changes transform the integral into the following:
genomic fragments. Environ. Microbiol. 6:938 –947.

(2␲)
1⫺M
2 冕 冘 e⫺ x2
i i dx
11. Krause, L., N. N. Diaz, A. Goesmann, S. Kelley, T. W. Nattkemper, F.
Rohwer, R. A. Edwards, and J. Stoye. 2008. Phylogenetic classification of
short environmental DNA fragments. Nucleic Acids Res. 36:2230 –2239.
12. McHardy, A., H. Martín, A. Tsirigos, P. Hugenholtz, and I. Rigoutsos.
BM⫺1[(L · K)1⁄2]
2007. Accurate phylogenetic classification of variable-length DNA frag-
where BM ⫺ 1[(L · K)1/2] is the M ⫺ 1 ball with a radius of (L
· K)1/2. ments. Nat. Methods 4:63–72.
After hyperspherical coordinate change, the integral becomes the 13. Blaisdell, B. 1986. A measure of the similarity of sets of sequences not
gamma distribution with (M ⫺ 1)/2 degrees of freedom requiring sequence alignment. Proc. Natl. Acad. Sci. U. S. A. 83:
L·K
5155–5159.

冊冕
14. Vinga, S., and J. Almeida. 2003. Alignment-free sequence comparison: a

冉M⫺1
冊 1 M⫺1
⫺s
review. Bioinformatics 19:513–523.
2 ,L·K ⫽
冉
p M⫺1 s 2 e ds 15. Karlin, S., and C. Burge. 1995. Dinucleotide relative abundance
⌫ 2
extremes: a genomic signature. Trends Genet. 11:283–290.
0
16. Sims, G., S.-R. Jun, G. Wu, and S.-H. Kim. 2009. Alignment-free
Thus, we let 1 ⫺ p[(M ⫺ 1)/2, L · K] be the P value for the com- genome comparison with feature frequency profiles (FFP) and optimal
parison of a query and a target genome under the hypothesis that they resolutions. Proc. Natl. Acad. Sci. U. S. A. 106:2677–2682.
17. Wu, T., Y. Hsieh, and L. Li. 2001. Statistical measures of DNA sequence
are related; i.e., a low P value indicates that the hypothesis that the dissimilarity under Markov chain models of base composition. Biometrics
query sequence originated from the target genome should be re- 57:441– 443.
jected. Notice that this P value is not the same as the Karlin-Altschul 18. Karlin, S., and S. Altschul. 1990. Methods for assessing the statistical
E value from sequence alignment, which measures the significance of significance of molecular sequence features by using general scoring
a score, obtained from a random distribution. schemes. Proc. Natl. Acad. Sci. U. S. A. 87:2264 –2268.
19. Greenbaum, B., A. Levine, G. Bhanot, and R. Rabadan. 2008. Patterns
ACKNOWLEDGMENTS of evolution and host gene mimicry in influenza and other RNA viruses.
PLoS Pathog. 4(6):e1000079.
This work has been supported by grants from by the Northeast Biodefense 20. Greenbaum, B., R. Rabadan, and A. Levine. 2009. Patterns of oligonu-
Center (grant U54-AI057158) and the National Library of Medicine cleotide sequences in viral and host cell RNA identify mediators of the host
(grant 1R01LM010140-01). innate immune system. PLoS One 4(6):e5969.

8 mbio.asm.org July/August 2010 Volume 1 Issue 3 e00156-10