Previously [20,21], we have described alignment free approaches to the problem of comparison and analysis of complete genomes, and some techniques enabling to cope with the sparseness of the n-gram type of genomic information representations. The problem of sparse occurrence matrices is even more pronounced when dealing with the number of permutations of the possible synonymous codons. Calculating the set of n–grams for such occurrences will lead to vector representations, which are severely sparse, especially for higher n–grams lengths, and hence to very poor statistics. To alleviate this problem a hybrid approach, [17] was proposed: namely, when computing counts of codon-pair patterns–separated by codon sub–sequences of differing length–the actual composition of these spacer sub–sequences will be neglected. However, when such partial counts are used as a composite set, poor statistics is no longer a hindering obstacle, and the complete information about particular n–gram frequencies profile is preserved, albeit in a distributed and convoluted form.

For every protein coding gene, with its original nucleotide sequence j₀, a set of equivalent nucleotide strings (j₁, j₂, j_3,…, j_N) is created by a Monte Carlo approach.

These artificial sequences have the following properties:

• • they are all of the same nucleotide lengths as the j₀;
• • they have exactly the same amino acid sequence as the j₀ (ie, the proteins translated from the j₁, j₂, j_3,…, j_N are identical to j₀);
• • they have, in the vast majority of cases, a synonymous codon order different from the original sequence j_0.

The last is an essential point, which merits a commentary. The probability that a given string j_i generated stochastically has the same synonymous codon order as the original j₀ decreases with the product of its length, with a probability limit tending rapidly to zero.

Therefore, the ISSCOR method allows comparing the original codon sequence with an ensemble of different synonymous sequences–yet all of them coding for the same sequence of amino acids.

Two sets of the H3N2 influenza virus hemagglutinin (HA) sequences were collected from NCBI: the set A comprising 2217 full length HA gene sequences (1701 nucleotides each); and the set B of 1810 sequences of the HA-1 antigenic region of the H3N2 hemagglutinin gene (987 nucleotides each). Of these 1810 sequences, the first 251 were identical with those described by Smith et al. [22] and isolated in the period of 1968–2002, and the remaining 1559 were obtained by collecting sequences with accession numbers given by Russell et al. [23], spanning the period of 2002–2006. The compositions of the sets A and B were almost entirely different as only two virus isolates were the same in both. The ISSCOR descriptors ([17], the equation 3 there), were calculated for these two sets of sequences, using codon spacer values from λ = 0 to λ = 16, and creating two matrices (MA and MB) of 2217 and 1810 rows respectively, and 2448 columns each.

It was observed that the G + C content in H3N2 hemagglutinin sequences of the set A decreased linearly over the period of 40 years, by approx. 5.93%–which corresponds to the changes in the relative ratios of all four nucleotides, although by different amounts for each. The linear regression parameters of the equation:

Results of principal component analysis (PCA) for the matrix MA (vide infra), showing 2217 data points, are depicted on Fig. 1A. There is a clear timeline trend of the sequences from the years 1968–1970, which are located at the upper left corner, going down towards the minimum at about the years 2001–2002, and then up again to reach the years 2007–2008 at the end of the right arm of the curve. Therefore, to examine this trend more synthetically, the yearly clusters were considered separately for each year between 1968 and 2008. For each yearly cluster their centers, as well as the corresponding standard deviations, were calculated, and are shown on Fig. 1B. The horizontal and vertical bars are proportional in sizes to the respective standard deviations of PC-1 and PC-2 values. Noteworthy, there are three regions of rather increased variability: first, around the year 1968, then the years 2000–2003, and finally the year 2008, with some less diverse periods in between as well. It was found earlier [25], based on extensive analyses of influenza HA sequences of the H1 and H3 subtypes, that the evolution of H3N2 hemagglutinin included long intervals of mostly neutral sequence evolution without noticeable antigenic change showing an excess of synonymous over nonsynonymous substitutions, punctuated by shorter intervals of rapid evolution during which newly dominant lineages quickly displaced previously coexisting ones [26,27]. The stasis intervals showed rather uniform distribution of replacements over the whole HA sequence's length, not favoring epitope regions.

The analogous PCA results for the matrix MB, showing 1810 data points, are depicted on Fig. 1C, marking also members of the antigenic clusters, described in Smith et al. [22]. As was the case for the set A, the yearly clusters of the set B were analyzed in the same manner, details are shown in Table 1 for all the years 1968–2006. Of special interest is the location of the A/Fujian/411/2002 sequence on this plot (the red diamond at about PC1 = 5 and PC2 = 0) as this is the oldest of the six strains constituting the Fu02 cluster, and yet all remaining sequences of this cluster (each isolated in 2003) have their PC2 values below it. Holmes et al. concluded that one clade of H3N2 viruses present at least since 2000 had provided the hemagglutinin gene for all the H3N2 viruses sampled after the 2002–2003 influenza season, and that a reassortment event was the likely progenitor of the antigenicaly variant influenza strains that caused the A/Fujian/411/2002 epidemic of the 2003–2004 season [28]. It is possible that a significant factor of such adaptation involved optimizing the functional compatibility of reassorting segments. Such a phenomenon might be a possible explanation why the lineage leading to the FU02 antigenic type did not dominate the viral population until a few years after its initial appearance, coincidental with an HA reassortment event [22,29]. It should be pointed out that the yearly clusters described here, and the antigenic clusters of Smith et al. do differ, as the latter are defined according to vaccination seasonal changes implemented in the past, thus they always contain sequences from more than a single year.

Antibodies against hemagglutinin provide protective immunity to influenza virus infection, and the HA is therefore the primary component of influenza vaccines, and as the antigenic structure of HA changes significantly over time, the vaccine has to be updated to ensure adequate efficacy against emerging viral variants. The WHO network of influenza centers routinely characterizes the antigenic properties of influenza viruses using inhibition assays, which combined with sequential data of variability in the antigenic HA-1 domain of the HA, are necessary to select strains for use in the seasonal vaccines [30], [31], [32]. Smith et al. [22] used antigenic data from 35 years of influenza surveillance between 1968 and 2003, with the resulting antigenic dataset consisting of a table of 79 post-infection ferret antisera by 273 viral isolates, and 4215 individual HA inhibition (HI) measurements, and then constructed an antigenic 2D map, to determine the antigenic evolution of influenza A/H3N2 virus, using projection [33] from the obtained high-dimensional antigenic data hyperspace.

The accuracy of the predictions has shown that their map might serve as a possible target of an attempt to describe antigenic relationships on a basis of the ISSCOR descriptors. Partial least squares regression (PLS-regression) is a technique used to find relationships between two data sets (X and Y), utilizing a latent variable (LV) approach to modeling the covariance possibly present in these two spaces. Rather than finding hyperplanes of maximum variance between the response and independent variables, as is the case in the PCA-regression, it finds a linear model by projecting the predicted variables together with the observables to a newly constructed space. In this way trying to uncover the multidimensional direction in the X space, that explains the maximum multidimensional variance direction in the Y space [34].

Therefore, the PLS-regression of the antigenic cluster centers’ ISSCOR descriptors, on the 2D map of the Smith's antigenic clusters, was performed, and the results are shown on Fig. 2A. The model utilizing six LVs was found to be optimal (RMS = 0.12), considering that the regression model obtained with only five LVs was not sufficient to achieve prediction errors small enough. Table 2 lists fifty of the most contributing ISSCOR variables used by each of the six LVs. At the bottom of each column there are values of explained variance for the ISSCOR descriptors [X-matrix], and the antigenic map's cluster centers [Y-matrix]. It is noteworthy that only two major LV would already suffice to explain 95% of variance in Y space but at the same time all six LV are necessary to explain the respective variance in X–the 5th and the 6th LV are both contributing almost equally strong.

The same PLS model was then used to project positions of H3N2 hemagglutinin strains isolated during each of the respective years (Fig. 2B). Of interest is a rather wide spread of the year's 2005 sequences, some of which are apparently reversing the general trend observed earlier (Fig. 2A), and continued subsequently by the majority of strains isolated in the other 2003–2006 yearly seasonal clusters. Russell et al. ([23], their supplementary materials) had identified and discussed several types of biases possibly influencing the composition of the set B of the HA-1 sequences, however, none of the biases they describe could convincingly explain the behavior of the yearly 2005 cluster observed here (in fact, dispersion in the year 2005 appears to be quite small–cf. Fig. 1B). On the other hand, there were also reports of surprising heterogeneity; strains isolated during 1998–1999 in Germany were characterized by a rather high variance: 75% of viruses examined during that period had completely different genome compositions [35], and similarly high local variability was found in human H3N2 influenza A viruses collected between 1999 and 2004 in New York State [28]. Such intermixing of antigenically distinct strains between different clusters was already present in the Smith's 2D map ([22]), where five strains (Hong Kong/3/95, AF008755; Hong Kong/32/95, AY661185; Hong Kong/38/95, AF008759; Victoria/75/95, AY661188; and Netherlands/91/96, AY661194) assigned to the Bejing 1992 (Be92) cluster were located instead within the Wuhan 1995 (Wu95) one, and vice versa–four strains assigned to the Wu95 (Madrid/G252/93, AY661151; Netherlands/372/93, AF008834; Netherlands/399/93, AY661167; England/7/94, Z46405) were found to belong to the Be92. The same nine sequences are also intermingled in our results from both the PCA (Fig. 1C), and the PLS regressions (Fig. 2A)–with the four strains from Wu95 mixing with the main Be92 cluster as in [22]; and the five strains from Be92 forming additional, small, intermediate cluster, located between the both parental clusters–together with additional eight strains assigned originally to Wu95 (Finland/339/95, AY661184; Netherlands/271/95, AY661191; Finland/338/95, AF368436; Finland/381/95, AY661196; Hong Kong/49/95, AY661187; Lyon/2279/95, AY661192; Geneva/AI9509/95, AY661182; and. Brazil/8/96, AF008733). Similar intermingling was also observed between the Bangkok 1979 (Bk79) and the Singapore 1987 (Si87) clusters (Fig. 2A), and a single strain from the Hong Kong 1968 cluster (Hk68) shifted to the England 1972 group (En72).

Significant vaccine-epidemic mismatches did occur in the past: during 1997–1998 season the Sydney-type viruses replaced the dominant Wuhan-type strains [30], and in 2003–2004, Fujian-like viruses replaced in turn Sydney strains [31]. The evolution of influenza HA has been the subject of intensive research for several decades [22], [36–44]. Many characteristics of HA's sequence evolution can be captured by standard Markov chain substitution models, which assign equal fitness to all accessible amino acids at a site. Kryazhimskiy et al. [45] have shown, however, that such models strongly underestimate the number of homoplasic amino acid substitutions during the course of HA's evolution. They found that such events preferentially occur at positively selected epitopic sites, which suggests that the evolution of the influenza's HA is strongly affected by the long-term site-specific preferences for individual amino acids. Modeling antigenic shift in influenza can help to predict the efficiency of vaccines. The virus is known to exhibit sporadic but sudden jumps in an antigenic distance, generating new antigenic clusters that mismatch vaccines in use, which makes prediction of antigenic properties of such novel strains a major challenge [26,27,43]. This is exacerbated by a scarcity of detailed antigenic data, which is a major problem, despite the very valuable effort of Ndifon et al. ([40], supplementary materials there) who initiated a public repository of such HI data. The current PLS model, which offers quite accurate correlation with the corresponding 2D antigenic map, was calibrated on only the eleven HA clusters’ centers. Should the antigenic data be made available for e.g. each of the all 251 sequences, and/or for more recent strains, improved calibration ought to yield much better and significantly more accurate description of antigenic relationships between A/H3N2 co-circulating viral variants.

The influenza A/H1N1 serotype was first found in the Spanish flu pandemic victims of 1918–1920, The origins of the 1918 strain are not clear, although it has been suggested that this virus had an avian origin [46,47]. The H1N1 subtype was subsequently a prevailing cause of seasonal influenza A outbreaks (together with subtypes B, and C) all along until 1957, when it was replaced by the A/H2N2 subtype (Asian flu), and later on in 1968 by the A/H3N2 (Hong Kong flu)–although inactivated A/H1N1 was included in antiviral vaccines even during periods when it was no longer a dominant influenza strain [48,49]. In March 2009 reemergence of the A/H1N1 was reported, first in Mexico and US, but the novel variant pandemic rapidly started to be a major source of infection, replacing A/H3N2 through the world. The lack of similarity between the novel variant virus and its nearest relatives indicated that its gene segments have been circulating undetected for an extended period [50]. Its low genetic diversity suggests that the introduction into humans was a single event, or multiple events of similar viruses.

Two sets comprising orthologous genes of the H1N1 influenza's virus hemagglutinin (HA) were collected from NCBI in January 2010; all sequences were unique and in cases when two (or several) sequences were found to be identical, only the earliest one was included in their respective sets:

• • the set C, comprising full length HA gene sequences of 1698 nucleotides (1698NN): 846 sequences of seasonal human A influenza H1N1 (H1N1s_C), 17 sequences of the novel, swine-like variant (H1N1v_C), and 20 sequences of H1N1 isolated from pigs (H1N1sw_C);
• • the set D, comprising full length HA gene sequences of 1701 nucleotides (1701NN): 173 sequences of seasonal human A influenza H1N1 (subset H1N1s_D), 727 sequences of the novel, swine-like variant (subset H1N1v_D), and 174 sequences of H1N1 isolated from pigs (subset H1N1sw_D).

The ISSCOR descriptors ([17], the equation 3), were calculated for these two sets of sequences, using Codon spacer values from λ = 0 to λ = 16, creating two matrices (MC and MD) of 1074, and 883 rows, respectively, and 2448 columns each. The PCA results, scatter plots of PC-1 vs. PC-2 are shown on Fig. 3A (for the set C) and Fig. 3B (for the set D).

In the context of analyzing A/H1N1 hemagglutinin sequences one of the most interesting questions is the possible origin of the novel serotype. Antigenically the H1N1 v viruses were homogeneous and similar to North American swine H1N1 viruses, but distinct from seasonal human H1N1. An early cluster analysis of the origins of the new influenza H1N1 virus of Rabadan et al. [52] indicated, that while the gene segments PB1, PB2, PA, HA, NP and NS are related to swine H1N2, and swine H3N2 influenza viruses isolated in North America, the segments NA and M are related to swine influenza viruses isolated in Eurasia. These findings corroborated slightly earlier results of Garten et al. [50] based on maximum likelihood phylogenetic analyses. Of the 1698 NN sequences isolated from swine hosts, which nevertheless all belong notably to the H1N1 human-isolated gene clusters: the D00837 (A/swine/Cambridge/39) was found earlier ([53], their Fig. 2), ([54], their Fig 3) to be located on the human evolutionary lineage. The same was true for the FJ638306 (A/swine/NC/00573/2005), and FJ638298 (A/swine/IL/00685/2005) ([50], their Fig. S2D HA). It is therefore safe to assume that also the FJ611898 (A/swine/Minnesota/07002083/2007) belongs to the same category of A/H1N1 human virus carried over to a swine (Fig. 3AFig. 3). All three sequences: FJ638306, FJ638298, and FJ611898 were found to form jointly one small clade on the neighbor joining tree, corroborating our findings from the ISSCOR analysis.

On the other hand, comparison between the two sets C and D reveals some notable differences. Of the 1701NN sequences isolated from swine hosts there are two: the EU004452 (A/swine/Henan/01/06), and EU004444 (A/swine/Tianjin/01/04) [50,51] both belonging to the large human seasonal H1N1 flu cluster at the top-left in the Fig. 3B. There are also several examples of reverse case–the 1701NN sequences isolated from humans, but grouped together in predominantly swine clusters. This can be observed for the S62154 (A/Alma Ata/1417/84) belonging to the 1930–1945 cluster of swine sequences–Gibbs et al. [54,55] have found S62154 among recombinant hemagglutinin genes, with the A/swine/Iowa/15/30 (EU139823) as its closest neighbor, and suggested them as all being direct descendents of 1918 Spanish flu viruses. Noteworthy, the EU139823 sequence is also closest to the S62154 on our plot. All of the 727 sequences of the novel, swine-like virus variant (subset H1N1v_D) clustered together in one large, compact set (green dots at the top-left in Fig. 3B). The next cluster to the right comprises of 16 swine sequences, and one human sequence (FJ986620, A/Ohio/01/2007), all of them isolated in the US between 2003 and 2007. During the search for possible predecessors of the recent pandemic of swine-origin novel A/H1N1 variant, the FJ986620 was found [56] to be one of a few sporadic cases of infection of humans with a triple-reassortant swine influenza virus prior to March 2009.

Another such human infection by swine reassortant virus DQ889689 (A/Iowa/CEID23/2005) [52,57], found in Iowa study, is among five human A/H1N1 sequences co-clustered in a large 1985–2009 assembly near the bottom of the Fig. 3B; four other human H1N1 isolates located in this cluster are: CY024925 (A/Ohio/3559/1988) [52], L24362 (A/MD/12/1991) [52], U53162 (A/WI/4754/1994) [58], and U53163 (A/WI/4755/1994) [58]. Another three human sequences are located in the next large cluster grouping genes found in the years between 1961 and 1986: CY026139 (A/Wisconsin/301/1976), CY021957 (A/New Jersey/1976) and CY044365 (A/New Jersey/8/1976) [59]. Additionally, Aras et al. [59] have described that two closest relatives to the novel H1N1 human virus were found either in the US (AF455680, A/swine/Indiana/P12439/2000), or in Eurasia (GQ229269, A/swine/Hong Kong/415/2004) [60], both of which are of the A/H1N2 subtype, their sequences were therefore included in the current analysis as well–on Fig. 3B they are marked within the same 2003–2007 group (magenta triangles) as already described for the FJ986620 sequence, and indeed one of them (AF455680) is least distant from the compact cluster of the A/H1N1 2009 pandemic sequences. G.J.D. Smith et al. [60] have estimated for the novel variant H1N1 hemagglutinin gene that the time of most recent common ancestors (TMRCA) of outbreak samples was around 28 August 2008 (between 1st April 2008, and 2nd January 2009). However, their estimate of duration of unsampled diversity was as long as about 10 years (8.41 to 11.02 yrs.).

The set E comprising of 613 orthologous genes of the H5N1 avian flu virus’ hemagglutinin (HA), of 1707 nucleotides each–isolated from avian hosts during the 1996–2009 period–was acquired from NCBI in January 2010. All sequences were unique–in cases when two (or several) sequences were identical, only the earliest one was included. The ISSCOR descriptors ([17] equation 3) were calculated for this set of sequences, using codon spacer values from λ = 0 to λ = 16, to create the matrix ME of 613 rows and 2448 columns.

The PCA scatter plot of PC1 vs. PC2 for the matrix ME, showing all 613 data points, is depicted on Fig. 4. For nine countries, with the most numerous populations of isolates, their plots are presented as separate panels, there are also separate panels showing remaining strains obtained from Africa, Asia and Europe. Already the analogous results for the A/H1N1 have shown some territorial clustering in several regions (Figs. 3A and B), however, in case of H5N1 here the groupings are remarkably well defined and in many instances also rather well separated between themselves. This is especially evident for strains isolated in several Far Eastern, and South-East Asian countries: China, Vietnam, Hong Kong, Indonesia, and Thailand (Fig. 4). On the other hand, isolates from Europe co-cluster together independently of their country of origin, which is also true for African and Middle-East countries–all of which form the large gathering at around PC1 = -3, and PC2 = 0 position. In roughly the same location there are also points for several Asian countries (all of which are different from the already mentioned five): Afghanistan, Pakistan, Bangladesh, India, Mongolia, Kazakhstan, and Japan. There are also two exceptions from Asia (well separated blue triangles in the panel “Asia” on Fig. 4)–one sequence from Japan (A/duck/Yokohama/aq10/2003, AB212280), and one from Myanmar (A/chicken/Pyigyitagon/204/2006, AB474081), both of them collocated together within the bottom cluster from China. The former H5N1 strain was isolated from duck meat processed for human consumption, imported to Japan from Shandong Province, China in 2003. That virus was antigenicaly different from other H5 viruses, including the Hong Kong H5N1 viruses isolated from humans in 1997 and 2003 [61]. The latter strain was isolated during the 1st outbreak of bird flu in Myanmar in March 2006, and it was subsequently found to belong to the clade-7 (the WHO H5N1 evolution nomenclature) of the highly pathogenic avian influenza (HPAI) [62]. There is also one, rather unexpected case of the isolate from a Belgian crested eagle (A/crested_eagle/Belgium/01/2004, DQ182483; the green square, shown here for clarity at the edge of the top cluster in the panel “Vietnam”, Fig. 4). On a closer examination this highly pathogenic H5N1 strain turned out to have been isolated from smuggled eagles, confiscated at the national airport, after an attempt to illegally enter them into Belgium from Thailand by a traveler [63].

Phylogeography of H5N1 viruses was actively researched recently in the regions of Southeast Asia and Far East which were repeatedly subjected to outbreaks of bird flu in poultry, especially in Vietnam [64,65], Thailand [66], and China [67–70], with conclusions rather similar to ours, concerning the yearly origins and seasonal spread of the A/H3N2 viruses [23]. Wallace et al. studied the geographic diffusion of H5N1 following the migration paths of the virus, by a way of a genetic phylogeography of H5N1's HA and NA sequences, and shown that the Chinese province of Guangdong was the source of multiple H5N1 strains spreading at both regional and international scales. In contrast, Indochina appeared to be a regional sink, at the same time demonstrating bidirectional dispersal among localities within the region [71]. They have shown further [72] that the virus migration was filtered out at some international borders, like between China and Vietnam or Thailand, but not at e.g. the one between China and Japan. Their results were recently reanalyzed by Hovmoeller et al. [73] who alleged that using a single tree, and a single optimization path, misestimates the frequency of transmission events, and moreover that the use of a single tree can fail to detect possible transmission events.

For each of the 44 countrywide clusters, their Euclidean centers were found–on the basis of the respective ISSCOR-deviate PC values (Figs. S1 and S2, in the supplementary material, show positions of each center, as well as standard deviations of each cluster in both PC1 and PC2 dimensions, cf. the figure legends for explanations)–and then the sequences closest to each center were established. The outlier sequence from Japan (A/duck/Yokohama/aq10/2003, AB212280, vide infra) was also included in the resulting set. In order to compare the results of this study to a standard phylogenetic technique, we have also calculated the corresponding Neighbour Joining trees. Comparison between both representations had shown a high level of similarities, and the relationships among different clades on the NJ tree corresponded reasonably well to those found through the ISSCOR-PCA method.