Reticulation at its best — an example from the oaks

One particular case where networks turn out to be a versatile tool is the study of low-level evolutionary patterns. This is especially so when we leave the comfort zone of well-sorted molecular markers, and use more than a single individual per species. Our recently published data set on (mostly Mediterranean) oaks, provides a nice example of this.

Why so few people study oaks at the intra-generic level

Oaks are notoriously difficult to study because they don't bother too much about species boundaries (which can be more or less obvious) and – at one point – decided to not sort their plastids at all (and full plastomes, as I once saw for myself first-hand, won't help). Hence, all reasonable phylogenetic reconstructions of oak evolution have been based on genetic data from the nucleome. However, this imposes a new problem — the sequenced nuclear gene regions allow the recognition of the major lineages (which recently have been formalized), but the closer one comes to the species level the more difficult it is to resolve anything at all.

Even the famous ITS region, which includes the weakly constrained internal transcribed spacer ITS1, and the structurally quite constrained ITS2, and have been frequently advocated as plant barcodes, turns out to be a two-edged sword. Relationships between the major intra-generic lineages is relatively clear, the ITS is pretty divergent down to the species level, but at the individual level, one faces a intra-genomic divergence that often outmatches inter-species differentiation.

In some groups, like the most speciose and most widespread white oaks (sect. Quercus), identical ITS variants exist from individuals / species separated today by thousands of kilometers of ocean or icy wasteland. One possible explanation is that oaks have very large population sizes, and they are wind-pollinated, so that they have a high capacity to permanently homogenize their genepools. Plastids, on the other hand, are only transmitted via the large fruit, the acorns, and the main animal vector for distributing acorns, the jaybirds, are sedentary birds. Their backup-vector, the squirrels are known to hoard a lot of acorns in a single place, but not for migrating globally (unless we assist them).

Nonetheless, we readily notice that the intra-individual differentiation patterns appear not to be entirely random, and so in our study we moved to another nuclear multi-copy spacer known to be more variable than the ITS1 and ITS2 (hence, largely ignored by molecular phylogeneticists) — the 5S intergenic spacer (5S IGS). It didn't help too much for solving the white oak puzzle (in western Eurasia), but did give us new insights into the two other western Eurasian sections: Ilex and Cerris.

The 'host-associate' framework

A cloned 5S-IGS (or ITS) sequence is not a good OTU, because we are usually not interested in a clone phylogeny (a mere sequence genealogy), but in the phylogenetic relationships between the individuals or species carrying the cloned sequence variants: the nuclear spacer population. Even networks struggle with such data, and my colleague Markus Göker came up with the idea to treat this in the form of hosts, the individuals, and associates, the cloned sequences found in the individual (Göker & Grimm, BMC Evol. Biol. 2008 — open access). There are several options to transfer the primary clone (associate) data into individual (host) data.

Options that we tested for transferring associate data into host data.
CM = character matrix, DM = distance matrix. CMhosts, independent used were morphological matrices. ENT — entropy, FRQ — frequence, CON — strict consensus, MOD — modal consensus, and SIZ — sample size, are character transformations implemented in Markus' g2cef, PBC and MIN are distance transformations implemented in pbc (these and other little helper programmes can be found here).

Using three cloned (ITS) datasets, we found that for these data the "Phylogenetic Bray-Curtis" (PBC — see the next figure) distance transformation outperforms the other tested options.

Computation of the "Phylogenetic Bray-Curtis" distance. It's a modification of the Bray-Curtis dissimilarity using the minimum distance for each covered row/column instead absence/presence. H1/H2 = hosts with different sets of associates (A1–A6)

Incidental but interesting insights

Whenever I come into contact with such data I advise the use of the PBC distance transformation as the basis for the main individual-level network, but also to run the MIN distance transformation: MIN will just calculate the minimum inter-clone distance between the clone samples of two individuals, and use this as the inter-individual distance.

Neighbour-net using the MIN transformation

The MIN network (above) is quite bushy for these data, because we naturally have many shared 5S-IGS variants among individuals of the same species, but occasionally also shared by individuals of different species. Nonetheless, it visualizes some basic differentiation patterns in the clone sample: compare e.g. the coherent cluster 3, the crenata-suber lineage (the 'Cork Oaks') — all individuals share a pair of very similar to identical 5S-IGS clones; and the divergent cluster 4, the 'Vallonea' oaks — all individuals have different sets of clones, but uniuqe 5S-IGS variants separating them from all other Cerris oaks (long proximal edge bundle).

Furthermore, we have potential F1 hybrids (morphologically intermediate) in our sample, and such hybrids, e.g. tj08, should have very low (to zero) MIN distances with members of their parental lineages.

However, the PBC network (below) is as beautiful as it gets — I really love this transformation, as it always comes up with something usable and interpretable.

Neighbor-net based on PBC-transformed inter-individual distances. See Simeone et al. (PeerJ PrePrints 2018 — open access pre-print) for a discussion.

However, this network was a last minute addition, because a happy little "accident" happened along the way, and the networks we were working with and looking at while drafting the paper where not PBC networks, as I thought.

It happened this way. Also implemented in Markus' little helper program are AVG, the average inter-clone distance, and MAX, the maximum inter-clone distance. AVG and MAX don't result in a proper distance matrix, because the diagonal will be the average or maximum distance between the clones of a single individual, and not all-zero as it should be (for MIN it's always zero). [We discussed a few options to modify AVG and MAX to ensure a zero diagonal, but couldn't devise something that makes sense.]

However, the SplitsTree program didn't bother about an all-zero diagonal, so the AVG and MAX transformed distance matrices will produce a Neighbor-net. So, what I assumed were PBC networks were in fact AVG networks.

Neighbor-net based on AVG-transformed inter-individual distances.

It took me quite long to recognize this "error" because, in contrast to the AVG (and MAX) networks I looked at when we did the 2008 paper, the one for the oaks made a lot of sense. Notably, the suspected F1 hybrids were perfectly resolved spanning up according boxes, and the species aggregates (clusters) did make sense regarding the general geographic setting, the history of the region under study, and their morphology.

Same graph as above, highlighting known or potential F1-hybrids spanning up according boxes.

For these data (with a minimum of four clones available per individual, individuals covering all species, and including the entire range of the section in western Eurasia), the AVG network better shows the potential F1 hybrids (or introgrades) than the (more methodologically sophisticated) PBC network. However, the latter makes more sense regarding speciation processes and the history of the group (because, the distance is a "phylogenetic" version of the well-known Bray-Curtis distances).

A "cactus-oak" fusion graph depicting nuclear and plastid differentiation (and evolution) in Quercus Group Cerris.

Take-home message

First, it's always good to delegate work you can do by heart to somebody new to it! This forces its propagation, which is important. More importantly, though, one has ones preferences and established analysis pipelines, and they may have become restricted in scope. I mainly used the -a (AVG), -i (MIN) and -x (MAX transformation) options in the little helper program to quickly summarize some of the primary differentiation data — for example, individuals have identical clones (MIN = 0), intra-individual divergence may be higher or not than inter-individual (MAX intra-individual > MIN inter-individual), and individuals may have strongly divergent clones (high MAX). AVG was computed and tabulated but never cherished by me. I always looked at the MIN transformed networks, since this provides a valid distance matrix, but then ignored them. But I never again tried to infer a Neighbor-net based on AVG or MAX transformations after our 2008 paper.

Second, Neighbor-nets are so quick to infer that there is no resource- or logic-related reason to not just run whatever distance one has on hand or can easily establish. Maybe even the biologically less-sound will reveal some interesting aspect (there are a lot of biological arguments that can be put forward for dismissing AVG distances in favour of PBC distances)

Paper (pre-print) and open data
Simeone MC, Cardoni S, Piredda R, Imperatori F, Avishai M, Grimm GW, Denk T. 2018. Comparative systematics and phylogeography of Quercus Section Cerris in western Eurasia: inferences from plastid and nuclear DNA variation. PeerJ Preprints 6: e26995v1.
Primary data and analysis files are included in the Online Supplemantary Archive: Simeone et al., PeerJ Preprints, doi: 10.7287/peerj.preprints.26995v1/supp-4. (See Readme.txt included in the topfolder of the archive.) 

Multimedia phylogeny?

Evolutionary concepts have often been transferred to other fields of study, or derived independently in them, especially in anthropology in the broadest sense, covering all cultural products of the human mind. This includes phylogenetic studies of languages, texts, tales, artifacts, and so on — you will find many examples of such studies in this blog. One of the more recent applications has been to what is sometimes called multimedia phylogeny — the research field that "studies the problem of discovering phylogenetic dependencies in digital media".

I have noted before that phylogenetics in the biological sense is an analogy when applied to other fields, because only in biology is genetic information physically transferred between generations — in the other fields, cultural information transfer is all in the minds of the people, not in their genes (see False analogies between anthropology and biology). This analogy often becomes problematic when applied to other fields, because the practical application of bioinformatics techniques separates the informatics from the bio, and the mathematical analyses focus on trying to implement the informatics without any biological justification.

A recent paper that discusses the application of bioinformatics to multimedia phylogeny exemplifies the potential problems:

Guilherme D Marmerola, Marina A Oikawa, Zanoni Dias, Siome Goldenstein, Anderson Rocha (2017) On the reconstruction of text phylogeny trees: evaluation and analysis of textual relationships. PLoS One 11(12): e0167822.

The authors described their background information thus:

Articles on news portals and collaborative platforms (such as Wikipedia), source code, posts on social networks, and even scientific publications or literary works, are some examples in which textual content can be subject to changes in an evolutionary process. In this scenario, given a set of near-duplicate documents, it is worthwhile to find which one is the original and the history of changes that created the whole set. Such functionality would have immediate applications on news tracking services, detection of plagiarism, textual criticism, and copyright enforcement, for instance.

However, this is not an easy task, as textual features pointing to the documents' evolutionary direction may not be evident and are often dataset dependent. Moreover, side information, such as time stamps, are neither always available nor reliable. In this paper, we propose a framework for reliably reconstructing text phylogeny trees, and seamlessly exploring new approaches on a wide range of scenarios of text reusage. We employ and evaluate distinct combinations of dissimilarity measures and reconstruction strategies within the proposed framework.

So, their solution to the separation of bio from informatics is to try a range of techniques, none of which are based on any particular model of how phylogenetic changes might occur in text documents. All of these methods involve distance-based tree-building.

The essential problem, as I see it, is that without a model of change there is no reliable way to separate phylogenetic information from any other type of information. For example, similarity can arise from many sources, only some of which provide information about phylogenetic history — phylogenetic similarity is a form of "special similarity". In biology, other sources of similarity are usually lumped together as chance similarities, such as convergence, parallelism, etc. Without this basic separation of phylogenetic and chance similarity, it does not matter how many distance measures you use, or how many tree-building methods you employ — if you can't separate phylogeny from chance then you are wasting your time constructing a hypothetical  evolutionary history.

The authors' only saving grace is their claim that: "In text phylogeny, unlike stemmatology [the analysis of hand-written rather than digital texts], the fundamental aim is to find the relationships among near-duplicate text documents through the analysis of their transformations over time." The expectation, then, is that the phylogenetic similarity of the texts will be high, which will thus reduce the possibility of chance similarities. Sadly, it will also reduce the probability that the similarities will contain any phylogenetic information at all — this is the classic short-branches-are-hard-to-reconstruct problem in phylogenetics.

For digital texts, the authors employ three distance measures: edit distance, normalized compression distance, and cosine similarity. None of these are model-based in any phylogenetic sense (although the first one is used in alignment programs such as Clustal) — I have discussed this in the post on Non-model distances in phylogenetics. Their tree-building methods include: parsimony, support vector machines (a machine-learning form of classification), and random forests (a decision-tree form of classification). Once again, none of these is model-based in terms of textual changes.

A final issue is the insistence on trees as the model of a phylogeny. In stemmatology, for example, a network is a more obvious phylogenetic model, because hand-written texts can be copied from multiple sources. Indeed, this distinction plays an important role in the first application of phylogenetics to stemmatology (see the post on An outline history of phylogenetic trees and networks). Perhaps this is not an issue for "near-duplicate text documents", but it does seem like an unnecessary restriction. Moreover, one of the empirical examples used in the paper actually has a network history, which therefore does not match the authors' reconstructed tree.

Non-model distances in phylogenetics

Bioinformatics is sometimes divorced from biology. This happens when the ideas involved are solely computational ones, and are not biologically motivated. For example, Otu et al. (2003), when proposing "a new sequence distance measure for phylogenetic tree construction" point out that "it is worth noting that our distance measures do not use any evolutionary model". They seem to see this a a Good Thing whereas biologists might not.

Pairwise distances are actually a good case in point, because there are many ways to measure the similarity / distance between any two objects, and most of them have nothing whatever to do with biology. One method that was popular about 10 years ago was using computerized compression.

This is based on the notion of complexity — all objects are complex to one degree or another, but the information they contain can often be reduced (or compressed) to a smaller amount without losing anything essential. That is, the original information can be exactly reproduced from the compressed version. This notion is applied to computer files, for example, in the idea of zipping a file — depending on the file's content it may be possible to zip the file to much a smaller size for storage or transmission.

This idea is straightforward to apply to similarity (Cilibrasi & Vitányi 2005). One simply has to compress the two objects separately as well as compress the combination of the two objects (ie. their union or concatenation). If the two objects are identical then the compression of their union will be the same size as the compression of each object individually (= a distance of 0), because all of the information in one of the objects is redundant. If they have nothing in common then the compression of their union will be the sum of the sizes of the compression of each object individually (= a distance of 1).

This idea has been applied in a number of areas where similarity / distance is used to group objects together (ie. clustering or classification), including computer viruses (Wehner 2008), music (Cilibrasi et al. 2004), languages (Benedetto et al. 2002a), and genomics (Chen et al. 2000, Li et al. 2001, Otu et al. 2003).

This approach has not always been received with equanimity. For example, the paper by Benedetto et al. (2002a) has received considerable commentary (Khmelev & Teahan 2003a; Benedetto et al. 2003; Khmelev & Teahan 2003b), (Goodman 2002; Benedetto et al. 2002b), (Wang 2009).

In biology, this idea has been mostly ignored. Chen et al. (2000) and Li et al. (2001) based their idea of information on Kolmogorov complexity, which is "the shortest program on a universal computer that outputs one of the objects when the input is the other object". This complexity cannot be measured directly, and so it is approximated by using a file compression program (eg. gzip). Clearly, the approximation can only be as good as the compression program and its algorithms.

Otu et al. (2003) based their idea of information on Lempel & Ziv complexity, which is "related to the number of steps required by a production process that builds either of the two objects". They claim an exact solution to the measurement of this complexity for gene sequences, and propose several slightly different distances based on this. Furthermore, the claim to:

show that the proposed distance measures, which are based on the relative complexity between sequences imply the evolutionary distance between organisms ... As the proposed approach does not depend on multiple alignments we test the validity of the approach in two ways: we use simulated data to show that the proposed distance measures can reasonably be represented by a tree. We also show the superiority of the proposed method on existing techniques using this simulated data. Secondly, we look at how well the results generated by the proposed method agree with existing phylogenies.

The use of generic compression similarity simply demonstrates that phylogeny does produce similarity to some extent. However, phylogeny is based on what we might call a form of "special similarity", which is similarity solely due to the possession of shared derived character states. Other components of similarity, such as the similarity produced by parallelisms, convergences and reversals, do not contribute to an estimate of a phylogeny. Indeed, they will often be positively misleading. A general concept of similarity is inadequate for reconstructing a phylogeny except under the simplest circumstances. This was empirically demonstrated, for example, in liguistics by the results of the Computer-Assisted Stemmatology Challenge, in which the compression method CompLearn came dead last in the primary ranking of the different phylogenetic methods (but did okay in the secondary ranking).

In short, only some aspects of complexity are relevant for phylogenetics, and only some of the information contained in a genome is relevant for measuring phylogenetic similarity. The focus on information as a whole ignores the "bio" in bioinformatics.

References

Benedetto D, Caglioti E, Loreto V (2002a) Language trees and zipping. Physical Review Letters 88: 048702. [arXiv:cond-mat/0108530, 2001]

Benedetto D, Caglioti E, Loreto V (2002b) On J. Goodman’s comment to "Language trees and zipping". arXiv:cond-mat/0203275

Benedetto D, Caglioti E, Loreto V (2003) A reply to the comment by Khmelev and Teahan. Physical Review Letters 90: 089804.

Chen X, Kwong S, Li M (2000) A compression algorithm for DNA sequences and its applications in genome comparison. Proceedings of the Fourth Annual International Conference on Computational Molecular Biology (RECOMB), pp. 107-117. ACM Press, Tokyo.

Cilibrasi R, Vitányi PMB (2005) Clustering by compression. IEEE Transactions on Information Theory 51: 1523-1545.

Cilibrasi R, Vitányi P, de Wolf R (2004) Algorithmic clustering of music based on string compression. Computer Music Journal 28(4): 49-67.

Goodman J (2002) Extended comment on "Language trees and zipping". arXiv:cond-mat/0202383

Khmelev DV, Teahan WJ (2003a) Comment on "Language trees and zipping". Physical Review Letters 90: 089803.

Khmelev DV, Teahan WJ (2003b) Comment on the reply of Benedetto et al. arXiv:cond-mat/0303261

Li M, Badger JH, Chen X, Kwong S, Kearney P, Zhang H (2001) An information-based sequence distance and its application to whole mitochondrial genome phylogeny. Bioinformatics 17: 149-154.

Otu HH, Sayood K (2003) A new sequence distance measure for phylogenetic tree construction. Bioinformatics 19: 2122-2130.

Wang X-L (2009) Comment on "Language trees and zipping". arXiv:0903.3669

Wehner S (2008) Analyzing worms and network traffic using compression. Journal of Computer Security 15: 303-320. (arXiv:cs/0504045v1, 2007)

Networks can outperform PCA ordinations in phylogenetic analysis

In a previous blog post I noted that there are potentially Distortions and artifacts in Principal Components Analysis analysis of genome data, in which the Principal Components Analysis (PCA) ordination axes are not independent of each other. This happens when the samples form one or more monotonic gradients, in which case higher-dimensional axes of the ordination plot are curved and twisted relative to the lower axes. The so-called Arch Effect is the best-known example of this problem. The danger is that the ordination diagrams will be interpreted as showing biological patterns when these artifactual patterns are not in the original datasets.

This mathematical artifact can be dealt with by: (i) mathematically correcting the distortion; (ii) using a Non-metric Multidimensional Scaling (NMDS) ordination instead; or (iii) using a phylogenetic network to display the multidimensional data, rather than using an ordination. Here, I illustrate this third possibility.

Network analysis of a 1-dimensional gradient

In the previous blog post I used the following small dataset as an example of a genetic dataset with a 1-dimensional pattern of relationships among the taxa. Here the taxa form a series of overlapping groups, defined by the pattern of Cs.

The diagonal dataset.

In the previous post I also illustrated the PCA ordination of this dataset, which shows the classic Arch Effect, rather than the expected linear arrangement of samples in the scatter plot. Below is a table showing whether analyses using two different types of data-display network can effectively display all of the taxon groups for datasets of this type but with varying numbers of taxa. All analyses were performed using the SplitsTree program.

As you can see, the NeighborNet analysis successfully displays all of the relevant bipartitions in all of the datasets. On the other hand, the related Neighbor-Joining tree only succeeds up to a taxon size of six. Six taxa can be represented by a simple linear tree, which can display three non-trivial bipartitions plus six trivial partitions, which is all that is needed when each taxon consists of five Cs. (Note that the NeighborNet and the Median network are identical to this tree, as well, for six taxa or less.)

The next two figures show the NeighborNet graphs for taxon sizes of 10 and 20, respectively, as examples.

NeighborNet of the diagonal dataset for n=10.

NeighborNet for the diagonal dataset for n=20.

The table also shows that the Median network analysis is only successful for small numbers of taxa. For greater numbers of taxa, some of the bipartitions are not displayed, as listed in the table. For example, for a taxon size of 10, the split {2,3,4,5,6}{1,7,8,9,10} is not displayed, as illustrated in the next figure.

Median network for the diagonal dataset for n=10.

For even larger numbers of taxa, the Median network superimposes some of the taxa, even though they have distinct sequences, and thus it cannot display all of the bipartitions, by definition. For example, for a taxon size of 20, taxon_9 and taxon_10 are superimposed, as shown in the next figure.

Median network for the diagonal dataset for n=20.

Network analysis of a 2-dimensional gradient

In the previous blog post I used the following dataset as an example of a genetic dataset with a 2-dimensional pattern of relationships among the taxa. Here the taxa form two series of overlapping groups, one defined by the pattern of Cs and the other by the pattern of Ts.

The grid dataset.

In the previous post I also illustrated the PCA ordination of this dataset, which shows a complicated 3-dimensional version of the 2-dimensional pattern. On the other hand, the Median network perfectly reproduces the original 2-dimensional grid arrangement of the taxa, as shown next, so that the network diagram is exactly as expected.

Median network for the grid dataset.

Sadly, the NeighborNet analysis does not fare quite so well, as it produces a much more complicated arrangement of the taxa. This arrangement does approximately represent the original grid-like pattern contained in the data, but it is presented in an unnecessarily convoluted way.

NeighborNet for the grid dataset.

Conclusion

It appears that data-display networks can be used to display multidimensional information, in the same manner as an ordination analysis. These networks try to display all of the relevant bipartitions in the data, rather than simply displaying a scatter plot of the samples.

However, there are considerable differences in the ability of the various network analyses to display the patterns inherent in the original data. In particular, NeighborNet seems to be very effective for displaying 1-dimensional patterns, while Median networks are more effective for 2-dimensional patterns.

Finally, it is worth noting that David Bryant has already thought of comparing NeighborNet with a Multi-dimensional Scaling (MDS) ordination, using an empirical dataset (2006, Radiation and network breaking in Polynesian linguistics. In: Forster P, Renfrew C (eds) Phylogenetic Methods and the Prehistory of Languages. McDonald Institute Press, University of Cambridge, pp 111-118). He found that they approximately agreed with each other.

Distortions and artifacts in Principal Components Analysis for analysis of genome data

It is perfectly clear from examples in the genomics literature that there is a problem occurring with the use of Principal Components Analysis (PCA) ordinations, and that it is apparently going undetected by users. This problem involves a spurious second axis in the output graph that is a curvilinear function of the first axis. The danger is that the ordination diagrams will be interpreted as showing biological patterns (see François et al. 2010; Arenas et al. 2013), when these patterns are not in the original datasets. I illustrate this problem here using some very simple datasets, and point to some examples of it in the literature. I then suggest some possible solutions, including the use of phylogenetic networks.

Introduction

Genomic data are multivariate, in the sense that we have multiple characters (the nucleotide sequences) for multiple organisms. "Ordination" is the term used when we arrange the multivariate data in a rational order, which is then usually displayed as a scatter plot, with the points representing the organisms. The spatial relationships of the points on the graph then indicate how similar the organisms are to each other (based on their genomes) — close points are more similar to each other than are distant points, as shown by the example in the first figure.

Salmela et al. (2011), Figure 2d.
Each dot represents the genome of one Swedish person,
coloured according to their geographical location within northern Sweden.

There are many mathematical techniques for ordination, although they fall into two main groups: (i) eigenanalysis methods, and (ii) non-metric multidimensional scaling methods. These methods have a long history in some parts of biology, for example in ecology (Gauch 1982; Jongman et al. 1987; Kent & Coker 1992; Leps & Smilauer; Legendre & Legendre 2012), although they are virtually unknown in other parts. (Note: the eigenanalysis ordinations are sometimes confusingly called metric multidimensional scaling ordinations.)

The eigenanalysis method called Principal Components Analysis (PCA) was introduced by Patterson et al. (2006) to the study of genomic data, having previously been introduced into genetics by Menozzi et al. (1978)Cavalli-Sforza and Edwards (1964), and it has become quite popular in the literature since then. The basic idea is to summarize the data by calculating the correlation coefficient between each pair of organisms, and then to display as much as possible of the common pattern of variation along the first ordered axis of the graph, and as much of the remainder along the second axis, and so on. The different axes are expected to show independent patterns of genetic relationship among the organisms, with the first axis showing the most important pattern, the second axis the second-most important pattern, and so on.

However, ecologists have long known that the first two components of a PCA (and related eigenanalysis techniques such as correspondence analysis) often show an arching pattern, which is a mathematical artifact of the eigenanalysis rather than a real biological pattern in the original data (Gauch et al. 1977). This is most apparent when the sampling units cover a long ecological gradient and those samples at each end of the gradient have few species in common (Minchin 1987; see also this blog post). This is called the Arch Effect (or the Horseshoe Effect, or the Guttman Effect), in which the second axis is curved and twisted relative to the first axis, and thus the axes are not independent of each other.

The problem illustrated

Basically, the distortion involves the data pattern being displayed using one more axis than is needed for that data pattern. That is, a 1-dimensional pattern is displayed using 2 dimensions, instead, and a 2-dimensional pattern is displayed using 3 dimensions, etc. I will first illustrate this with a simple 1-dimensional pattern.

This first dataset consists of 20 nucleotide sequences. Each sequence consists of only two nucleotides (19 A and 5 C), and the pattern of Cs forms a simple gradient (ie. there is only one pattern in the data). This gradient can be displayed using only one ordination axis, which would arrange the taxa in order from Taxon_1 to Taxon_20. However, note that the patterns at the two ends of the gradient have nothing in common as far as the pattern of Cs is concerned.

The gradient dataset.

The next figure is the Principal Components Analysis (PCA) ordination diagram of these data. Note that the 1-dimensional pattern is displayed using 2 dimensions, arranging the 20 taxa in an arch rather than in a straight line. (Note that the open part of the arch can face either up or down; in this case it happens to face up.) So, it appears in the scatter plot that Taxon_1 and Taxon_20 are rather similar to each other (ie. they are near each other), which is not true. Furthermore, the dots are not equally spaced, which they should be given the pattern in the original data. The second dimension, along with the unequal spacing of the points, is a mathematical artifact of the eignenanalysis rather than a pattern in the data. It is caused by the unimodal distribution of the pattern of Cs along the gradient.

PCA ordination of the gradient data.
Not all of the taxa are labelled.

The second dataset consists of 24 nucleotide sequences. Each sequence consists of three nucleotides (8 A, 5 C and 4 T), and the pattern of Cs forms one gradient while the pattern of Ts forms a second gradient (orthogonal to the first one). These two independent patterns can be displayed using two ordination axes, which would arrange the taxa in a 6x4 rectangular grid. (If this is not clear, look at the bottom two figures of this post.)

The grid dataset.

This next figure is the PCA ordination diagram of these data. Note that the 2-dimensional pattern is displayed using 3 dimensions, only two of which are actually shown in the diagram. (It is standard to show only the first two axes, on the assumption that these contain most of the important information.) These axes arrange the 24 taxa in an arch in the first two dimensions, rather than in a straight line, and then separate pairs of taxa along the third axis. So, in two dimensions pairs of taxa appear to sit on top of one another, indicating a similarity that is not present in the original dataset. Once again, the use of the third dimension is an artifact of the eignenanalysis rather than a pattern in the data. It is caused by the unimodal distribution of the patterns of Cs and Ts along the two gradients.

PCA ordination of the grid data.

Clearly, when the samples form gradients then eigenanalysis is not an appropriate ordination technique, although it is very useful under many other circumstances. This problem with gradients has previously been pointed out for genomic data when they form a clinal geographical pattern (Novembre & Stephens 2008; Reich et al. 2008), but clearly the problem involves much more than just geographical patterns. Also, the problem is being ignored by users — most of the examples below were published after the paper by Novembre & Stephens.

Literature examples

The Arch Effect has appeared repeatedly in the genomics literature (often produced using the SmartPCA, Plink or Eigenstrat computer programs), and some examples are illustrated here. Note that the axes are not always in the same orientation, but the arch should be found on axis 2 relative to axis 1, and it can be oriented either up or down.

These first two examples come from the very paper that first introduced PCA to genomics studies.

Patterson et al. (2006), Figure 5.

Patterson et al. (2006), Figure 8.

These next two examples are from studies of genomic data from Scandinavia.

Salmela et al. (2011), Figure 2c.

Skoglund et al. (2012), Figure 1c.
Note that axis 1 is shown vertically.

The next two examples also come from worldwide genome studies. Note that the second of these refers to the ordination as metric multidimensional scaling (MDS).

Stoneking & Krause (2011), Box 1a.

Jakobsson et al. (2008), Figure 1d.

The next two examples also come from a worldwide genome study, but in this case the relationship of the Arch Effect to the geographical sampling is made explicit.

Wang et al. (2012), Figure 1b.

Wang et al. (2012), Figure 5b.

The next example illustrates rotation of the ordination axes, so that the arch is now aligned along the two axes.

Novembre and Ramachandran (2011), Figure 4e.

Finally, the Gene Expression blog has an example, where the ordination is called a metric multidimensional scaling (MDS) ordination. This is based on ~3,000 individuals and 250,000 SNPs.

Possible solutions

There are three possible solutions to the Arch Effect problem in eigenanalysis ordinations.

First, one can try to mathematically correct the distortion. This approach was actively pursued in ecology via what is known as detrending and rescaling (Hill & Gauch 1980), in which the axes are forced to be independent of each other after the eigenanalysis. This approach works, in the sense that the Arch Effect disappears, but it has the drawback that all relationships between the axes are removed, including those that require more than one dimension to be displayed (Palmer 1993). So, while this method is still commonly encountered in the ecological literature, it is not uniformly recommended.

This approach is illustrated in the next figure, which is a Detrended Correspondence Analysis (DCA, with detrending by second-order polynomials) ordination analysis of the second of my simple datasets above. Note that the original 2-dimensional grid arrangement of the taxa is perfectly preserved, so that the ordination diagram is exactly as expected.

DCA ordination of the grid data.

The second, and most obvious, solution to the problem is not to use an eigenanalysis ordination in the first place. That is, not all ordination techniques produce the same result. In this case we could use a Non-metric Multidimensional Scaling (NMDS) ordination, instead. The latter is a very different type of technique, mathematically, which explicitly tries to arrange the points in the ordination plot so that their relative distances reflect as closely as possible the pairwise "distances" in the original dataset (here, a distance is the inverse of a similarity). This technique is known to avoid the Arch Effect (Minchin 1987).

This is illustrated in the next figure, which is a NMDS ordination analysis of the second dataset. Note that the original 2-dimensional grid arrangement of the taxa is almost perfectly preserved. The main disadvantages of NMDS are the computational time required, the difficulty of finding the optimal solution for large datasets (eigenanalysis is significantly faster; Zheng et al. 2012), and the fact that the number of ordination axes needs to be specified beforehand.

NMDS ordination of the grid data.
Not all of the taxa are labelled.

The third solution to the problem is to use a phylogenetic network to display the multidimensional data, rather than using an ordination. Conceptually, this is similar to the NMDS type of ordination, but it uses a line graph to display the results rather than using a scatter plot. I have illustrated this in the next blog post (Networks can outperform PCA ordinations in phylogenetic analysis), where I report that this solution works quite well, too.

Conclusion

Every time you see an arch in a PCA plot, you might very well be looking at a mathematical artifact. You should bear this in mind when you interpret the ordination pattern, especially if you wish to draw biological conclusions from the diagram. You will, of course, see plenty of PCA ordinations without any evidence of an arch, and these are probably safe to interpret naively.

PCA can be useful when we expect genomes to be monotonically related to one another, but this is likely to be rather rare in nature. It is much more likely that genomes will have nonlinear relationships, such as unimodal patterns, with a peak at an intermediate part of any gradients that are present. Under these circumstances it produces artifacts, which may be mistakenly interpreted as representing real biological patterns. Under these circumstances, an alternative method of multivariate analysis should be used.

More information about ordination techniques can be found at the Ordination Methods web page.

Note: updated Friday December 21 2012.

References

Arenas M, François O, Currat M, Ray N, Excoffier L (2013) Influence of admixture and paleolithic range contractions on current European diversity gradients. Molecular Biology and Evolution 30: 57-61.

Cavalli-Sforza L, Edwards AWF (1964) Analysis of human evolution. Proceedings of the 11th International Congress of Genetics, The Hague, 1963. Genetics Today 3: 923-933.

François O, Currat M, Ray N, Han E, Excoffier L, Novembre J (2010) Principal component analysis under population genetic models of range expansion and admixture. Molecular Biology and Evolution 27: 1257-1268.

Hill MO, Gauch HG (1980) Detrended Correspondence Analysis: an improved ordination technique. Vegetatio 42: 47-58.

Gauch HG (1982) Multivariate Analysis in Community Structure. Cambridge University Press, Cambridge.

Gauch HG, Whittaker RH, Wentworth TR (1977) A comparative study of reciprocal averaging and other ordination techniques. Journal of Ecology 65: 157-174.

Jongman RHG, ter Braak CJF, van Tongeren OFR (1987) Data Analysis in Community and Landscape Ecology. Pudoc, Wageningen.

Kent M, Coker P (1992) Vegetation Description and Analysis: a Practical Approach. Belhaven Press, London.

Legendre P, Legendre L (2012) Numerical Ecology, 3rd edition. Elsevier, Amsterdam.

Leps J, Smilauer P (2003) Multivariate Analysis of Ecological Data using Canoco. Cambridge University Press, Cambridge.

Menozzi P, Piazza A, Cavalli-Sforza L (1978) Synthetic maps of human gene frequencies in Europeans. Science 201: 786-792.

Minchin PR (1987) An evaluation of relative robustness of techniques for ecological ordinations. Vegetatio 69: 89-107.

Novembre J, Ramachandran S (2011) Perspectives on human population structure at the cusp of the sequencing era. Annual Review of Genomics and Human Genetics 12: 245-74.

Novembre J, Stephens M (2008) Interpreting principal component analyses of spatial population genetic variation. Nature Genetics 40: 646-649.

Palmer MW (1993) Putting things in even better order: the advantages of Canonical Correspondence Analysis. Ecology 74: 2215-2230.

Patterson N, Price AL, Reich D (2006) Population structure and eigenanalysis. PLoS Genetics 2: e190.

Reich D, Price AL, Patterson N (2008) Principal component analysis of genetic data. Nature Genetics 40: 491-492.

Salmela E, Lappalainen T, Liu J, Sistonen P, Andersen PM, Schreiber S, Savontaus ML, Czene K, Lahermo P, Hall P, Kere J (2011) Swedish population substructure revealed by genome-wide single nucleotide polymorphism data. PLoS One 6: e16747.

Skoglund P, Malmström H, Raghavan M, Storå J, Hall P, Willerslev E, Gilbert MT, Götherström A, Jakobsson M (2012) Origins and genetic legacy of Neolithic farmers and hunter-gatherers in Europe. Science 336: 466-469.

Stoneking M, Krause J (2011) Learning about human population history from ancient and modern genomes. Nature Reviews Genetics 12: 603-614.

Wang C, Zöllner S, Rosenberg NA (2012) A quantitative comparison of the similarity between genes and geography in worldwide human populations. PLoS Genetics 8: e1002886.

Zheng X, Levine D, Shen J, Gogarten SM, Laurie C, Weir BS (2012) A high-performance computing toolset for relatedness and Principal Component Analysis of SNP data. Bioinformatics 28: 3326-3328.

Eurovision Song Contest 2006: a network analysis

Data-display networks can be used for displaying affinities between any group of objects, especially those for which a distance matrix can be calculated. Their use is essentially as Exploratory Data Analysis (EDA), creating a visual display of the patterns of relationship (or affinity), without strong assumptions about the cause of the affinities.

Philippe Gambette has, tucked away on a web page linked from his blog, a brief tree-based analysis of the scores for the finalists in the Eurovision Song Contest from 2006. It seems to me that these data illustrate some interesting things about the possible uses of networks in EDA, and so I will expand on the analysis here.

For those of you who know little of this odd musical institution (held every year for a bit longer than I've been alive), I will provide a brief description. All countries in Europe (very broadly defined, being any country that is a member of the European Broadcasting Union) are entitled to submit a song for consideration, although some countries have never participated (and others have rarely contributed). The submitted songs are performed, and are voted upon by all of the countries participating that year.

There were 24 countries represented in the Final in 2006: Armenia, Bosnia & Herzegovina, Croatia, Denmark, Finland, France, Germany, Greece, Ireland, Israel, Latvia, Lithuania, the Former Yugoslav Republic of Macedonia, Malta, Moldova, Norway, Romania, the Russian Federation, Spain, Sweden, Switzerland, Turkey, the Ukraine and the United Kingdom. The entries from Albania, Andorra, Belarus, Belgium, Bulgaria, Cyprus, Estonia, Iceland, Monaco, the Netherlands, Poland, Portugal and Slovenia did not make it into the final, but these 13 countries retained their voting rights. Serbia & Montenegro withdrew from the contest but also retained voting rights, resulting in 38 voting countries.

Red = did not qualify for the final; green = qualified for the final; orange = non-participating in 2006; grey = never participating.

Voting was carried out in 2006 by having each country assign scores of 1, 2, 3, 4, 5, 6, 7, 8, 10 and 12 to ten of the 24 contestants in the Final. The winner was the one who accumulated the greatest total. As it turned out, these totals ranged from 1 to 292 per song that year.

Network analysis

Philippe Gambette suggests the Manhattan distance as a suitable measure of the relationships in voting patterns among the 38 voting countries. I will discuss this choice below, but for the moment let's proceed with this choice. I have visualized the distance matrix with a NeighborNet analysis. On this network I have superimposed some broad geographical groups as different colours. Countries that are closely connected in the network are similar to each other based on their voting patterns, and those that are further apart are progressively more different from each other.

NeighborNet based on the Manhattan distance. Click to enlarge.

I am sure that you can recognize the geographical groupings. Red represents the countries from Northern Europe and around the Baltic; green is for the countries from Eastern Europe; orange is for the countries of Western Europe; blue is for the countries from Southern Europe and the Middle East; and purple is for the countries from the former Soviet Union, in Far Eastern Europe. Clearly, given these definitions, some of the larger countries could be in more than one group (e.g. France, Poland); and I could, of course, have chosen any way of grouping the countries that I liked.

However, I feel that it is noteworthy that the voting follows geographical lines, with few exceptions. Obviously, it is not geography per se that is reflected in the voting, but sociological influences associated with the geographical distributions of the cultural, political and language groups (see Clerides & Stengos 2006; Raykoff & Tobin 2007; Ginsburgh & Noury 2008; Spierdijk & Vellekoop 2009).

This network is EDA and thus only a data summary, and so we should look at the patterns in more detail, to see how they are created.

The Western European countries (as defined in my groupings) gave their big votes to Turkey, Finland, Bosnia & Herzegovina, Armenia and Greece (in that order). The Eastern European countries gave their big votes to Bosnia & Herzegovina, Russia, Finland, Greece, Macedonia and Croatia. Three of the countries are held in common between these two lists, with Finland highest.

The Southern Europeans gave their big votes to Romania, Russia, Finland and the Ukraine. The former Soviet countries gave their votes to Russia, Bosnia & Herzegovina, the Ukraine, Sweden and Lithuania. Here, two countries are in common between these two lists, with Russia highest.

The Northern European countries gave their big votes to Finland, Lithuania, Russia and Sweden, in that order.

Note that the Finns (who won) appear on four of these five lists, while Russia (who came second) also appears on four of the five. However, Finland does not appear on the former Soviets' list, and Russia does not appear on the Western European list.

This also makes it clear why the Finnish vote does not follow geographical lines: Denmark, Estonia, Greece, Iceland, Norway, Poland, Sweden and the United Kingdom, for example, all gave their 12-point vote to Finland, whereas Finland (who couldn't vote for themselves) gave their main votes to Russia and then Bosnia & Herzegovina. It is thus important to note that the winner's voting pattern will rarely be "geographical", because the competition regulations over-rule this possibility.

Andorra (an outlier in 2006) gave their big points to Spain (which no-one else did) but the rest of their voting exactly followed the Northern European countries. The only other unique vote was the 12 points that Malta gave to the Swiss song.

Armenia straddles two geographical groups in the network because they gave their big points to Russia and the Ukraine, but most of their other points went to Eastern European countries. Monaco has a long terminal edge because they gave many votes to Ireland and Latvia, which few other countries did.

Austria, Georgia, the Czech Republic and Hungary all announced that they would not be participating in 2006, and Italy had not taken part in the Contest since 1997 (some of these countries are orange on the map above). Other qualified countries had never previously participated (and some still have not) (grey on the map). Based on the EDA analysis, however, we might be able to predict what their voting pattern would be.

For example, we might predict that the Czech Republic, Hungary and Slovakia would vote with the Eastern European countries, while Azerbaijan and Georgia would vote with their former Soviet compatriots. Austria, Liechtenstein and Luxembourg are likely to have voted with the Western Europeans, while Italy and San Marino (and the Vatican City!) would presumably be more influenced by the Southern European choices.

Alternative analysis

As noted above, an interesting question for EDA is the choice of distance measure. In particular, we need to consider whether double-zeroes (negative matches) should count as representing a similar voting pattern between two countries. Given that each country can vote for only 10 songs out of the 24 available, there will be between 4 (if there are no votes in common) and 14 (if they vote for the same ten songs) countries that are not voted on by any given pair of countries. Do these "double absent" votes count as a similar voting pattern or not? The above analysis assumes that they do, but this may not be a reasonable assumption, because the joint absence of at least 4/24 votes is a requirement of the rules rather than a natural choice of the voters. Once again, the competition regulations over-rule the possibility of geographical patterns — we simply do not find out how each of the countries would have voted beyond 10/24 finalists (or 10/23 when the voter is also a competitor).

This is a clear case where absence of evidence is being treated as evidence of absence. This confounding of two sources of "experimental" variation seems to be a factor in most of the previous analyses of the Eurovision voting patterns that I have been able to locate, including those based on clustering (eg. Yair 1995; Fenn et al. 2006; Dekker 2008), ordination (eg. Yair 1995; Yair & Maman 1996; Doosje & Haslam 2005; Dekker 2008), multiple linear regression (eg. Bruine de Bruin 2005; Haan et al. 2005; Ginsburgh & Noury 2008; Spierdijk & Vellekoop 2009), z-scores (eg. Gatherer 2004), decision trees (Ochoa et al. 2009) and monte carlo simulations (eg. Gatherer 2006). In all cases, lack of a vote is treated as a score of zero rather than treated as "not applicable", so that absence of evidence (not applicable) is being treated as evidence of absence (a zero score). The only exception I found is the use of a censored regression model, such as the rank-ordered logit or tobit used by Clerides & Stengos (2006), where "not applicable" is treated as a censored observation.

An alternative distance measure that ignores all double-zeroes is the Steinhaus dissimilarity (or Bray-Curtis similarity), which is basically the Manhattan distance standardized by the total of the scores for each pair of countries. Using this distance measure means that at least one of each pair of countries must have voted for a song before it counts as a measurement of the voting patterns. This does not deal with the issue of comparing a "not applicable" result with a vote, but does deal with the issue of comparing two "not applicables" with each other.

NeighborNet based on the Steinhaus distance. Click to enlarge.

Using this distance instead does, indeed, have a significant effect on the observed patterns of voting relationships shown by the NeighborNet analysis. Much of the geographical pattern is still evident in the network, but it is fragmented into subgroups. Note, in particular, that Poland (with its ambiguous geographical placement) changes its affinities but France does not; and most importantly, Finland nows joins its geographical neighbours, because the disruptive effects of the rules have been eliminated from the analysis.

It is likely that the non-geographical patterns are related to some of the other factors that have been implicated to influence the Eurovision scores, such as performance order (Bruine de Bruin 2005; Haan et al. 2005; Clerides & Stengos 2006), host country (Bruine de Bruin 2005; Clerides & Stengos 2006), song language (Clerides & Stengos 2006) and possibly even song quality (Clerides & Stengos 2006; Ginsburgh & Noury 2008).

The network analysis is thus an effective EDA tool because it has revealed confounded experimental patterns that have not been addressed by the previous statistical analyses.

It is important to note that the difference between the Manhattan and Steinhaus distances is that the former considers patterns based on which countries are voted for as well as which ones are not, whereas the latter considers only the positive voting patterns. Since these distances produce different EDA data summaries, we can conclude that the strong "geographical" pattern shown by the Manhattan distance is, to a large extent, controlled by which countries are avoided in the voting rather than which are chosen to receive scores. I have long suspected that European politics works in precisely this same manner!

References

Bruine de Bruin W. (2005) Save the last dance for me: unwanted serial position effects in jury evaluations. Acta Psychologica 118: 245-260.

Clerides S., Stengos T. (2006) Love thy neighbor, love thy kin: Voting biases in the Eurovision Song Contest. Department of Economics, University of Cyprus, Discussion Paper 2006-01.

Dekker A. (2008) The Eurovision Song Contest as a ‘friendship’ network. Connections 28(1): 59-72.

Doosje B., Haslam S.A. (2005) What have they done for us lately? The dynamics of reciprocity in intergroup contexts. Journal of Applied Social Psychology 35: 508-535.

Fenn D., Suleman O., Efstathiou J., Johnson N.F. (2006) How does Europe make its mind up? Connections, cliques, and compatibility between countries in the Eurovision Song Contest. Physica A: Statistical Mechanics and its Applications 360: 576-598.

Gatherer D. (2004) Birth of a meme: the origin and evolution of collusive voting patterns in the Eurovision Song Contest. Journal of Memetics - Evolutionary Models of Information Transmission 8: 4.

Gatherer D. (2006) Comparison of Eurovision Song Contest simulation with actual results reveals shifting patterns of collusive voting alliances. Journal of Artificial Societies and Social Simulation 9: 2.

Ginsburgh V., Noury A.G. (2008) The Eurovision Song Contest: Is voting political or cultural? European Journal of Political Economy 24: 41-52.

Haan M., Dijkstra G., Dijkstra P. (2005) Expert judgment versus public opinion: Evidence from the Eurovision Song Contest. Journal of Cultural Economics 29: 59-78.

Ochoa A., Muñoz-Zavala A.E., Hernández-Aguirre, A. (2009) A hybrid system approach to determine the ranking of a debutant country in Eurovision. Journal of Computers 4: 713-720.

Raykoff I., Tobin R.D. (eds) (2007) A Song for Europe: Popular Music and Politics in the Eurovision Song Contest. Ashgate, Farnham UK.

Spierdijk L., Vellekoop M. (2009) The structure of bias in peer voting systems: Lessons from the Eurovision Song Contest. Empirical Economics 36: 403-425.

Yair G. (1995) 'Unite Unite Europe': The political and cultural structures of Europe as reflected in the Eurovision Song Contest. Social Networks 17(2): 147-161.

Yair G., Maman D. (1996) The persistent structure of hegemony in the Eurovision Song Contest. Acta Sociologica 39: 309-325.