Phylogenetic Trees from Large Datasets

Inaugural – Dissertation
zur
Erlangung des Doktorgrades derMathematisch–Naturwissenschaftlichen Fakultätder Heinrich–Heine–Universität Düsseldorf

vorgelegt von
Heiko A. Schmidt
aus BonnDüsseldorf
2003

Gedruckt mit der Genehmigung der Mathematisch-Naturwissenschaftlichen
Fakultät der Heinrich-Heine-Universität Düsseldorf

Referent:
Prof. Dr. Arndt von Haeseler
Korreferent: Prof. Dr. William Martin
Tag der mündlichen Prüfung: 28.07.2003

The Dachshund
Indeed, it was made on the plan of a bench for length and lowness. It seemed
to be satisfied, but I thought the plan poor, and structurally weak, on accountof the distance between the forward supports and those abaft. With age the
dog’s back was likely to sag; and it seemed to me that it would have been astronger and more practicable dog if it had had some more legs.
Mark Twain (Following the Equator, 1895-1896)

Biologists must constantly keep in mind that what they see was not designed,
but rather evolved.
Francis Crick (What Mad Pursuit, 1988)

iv

PrefaceThis thesis deals with topics from Bioinformatics/Phyloinformatics. It touches the fields
algorithmic complexity, molecular phylogenetics, sequence evolution, as well as parallel
computing and, thus, is addressed to an interdisciplinary community. Although it ishardly possible to write a completely self-contained thesis, I tried to provide some basic
information needed for members of either community to be able to more easily understand the topics discussed. Consequently, some information might be enclosed that seemtrivial to, e.g., computer scientists but which biologists are possibly not familiar withand vice versa.
When writing up my thesis I decided, as it is custom in scientific presentations, to use
the scientific ’we’, which I prefer over the personal pronoun ’I’.
Parts of this thesis have been published in the following articles:
1. H. A. Schmidt, K. Strimmer, M. Vingron, and A. von Haeseler (2002) TREEPUZZLE: Maximum likelihood phylogenetic analysis using quartets and parallel
computing. Bioinformatics, 18, 502-504.
2. H. A. Schmidt and A. von Haeseler (2002) Quartet Trees as a Tool to Reconstruct
Large Trees from Sequences. In K. Jajuga, A. Sokolowski, and H.-H. Bock (eds.),
Data Analysis, Classification, and Related Methods, Studies in Classification, Data
Analysis, and Knowledge Organization, pages 379-388, Springer, Heidelberg, NewYork.
3. H. A. Schmidt and A. von Haeseler (2003) Maximum Likelihood Analysis Using
TREE-PUZZLE. In A. D. Baxevanis, D. B. Davison, R. D. M. Page, G. Stormo,
and L. Stein (eds.), Current Protocols in Bioinformatics, Unit 6.6, pages 6.6.16.6.25, Wiley and Sons, New York.
4. H. A. Schmidt, E. Petzold, M. Vingron, and A. von Haeseler (2003) MolecularPhylogenetics: Parallelized Parameter Estimation and Quartet Puzzling. J. Parallel Distrib. Comput., 63, in press.
The TREE-PUZZLE package including developments presented in this thesis is freely
available from http://www.tree-puzzle.de.

Contents

1. Overview
1.1. Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.2. Organization of this Thesis . . . . . . . . . . . . . . . . . . . . . . . . . .

1
1
2

2. Molecular Phylogenetics: A General Introduction
2.1. Biological Sequences and Molecular Evolution . . . . . . . . . . . . . . .
2.1.1. Types of Biological Data . . . . . . . . . . . . . . . . . . . . . . .
2.1.2. Sequence Evolution . . . . . . . . . . . . . . . . . . . . . . . . . .
2.1.3. Multiple Sequence Alignment . . . . . . . . . . . . . . . . . . . .
2.1.4. Modeling Molecular Evolution . . . . . . . . . . . . . . . . . . . .
2.2. Phylogenetic Tree Reconstruction . . . . . . . . . . . . . . . . . . . . . .
2.2.1. Notation on Phylogenetic Trees . . . . . . . . . . . . . . . . . . .
2.2.2. Types of Phylogenetic Methods . . . . . . . . . . . . . . . . . . .
2.2.3. Maximum Likelihood on Phylogenetic Trees . . . . . . . . . . . .
2.2.4. Complexity of Phylogenetic Analysis . . . . . . . . . . . . . . . .

5
5
5
6
7
7
10
10
14
14
17

3. The Quartet Puzzling Algorithm
193.1. Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
3.2. Methods from the TREE-PUZZLE Package . . . . . . . . . . . . . . . . 20
3.2.1. Likelihood Mapping . . . . . . . . . . . . . . . . . . . . . . . . . . 20
3.2.2. Quartet Puzzling . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
3.2.2.1. ML Step . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
3.2.2.2. Puzzling Step . . . . . . . . . . . . . . . . . . . . . . . . 23
3.2.2.3. Consensus Step . . . . . . . . . . . . . . . . . . . . . . . 24
4. Improvement in Complexity of the Puzzling Step of QP
25
4.1. Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
4.2. Complexity Measures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
4.3. Complexity of the Original Puzzling Step . . . . . . . . . . . . . . . . . . 27
4.4. More Efficient Puzzling Step Algorithms . . . . . . . . . . . . . . . . . . 29
4.4.1. Split-Based Puzzling Step Algorithm . . . . . . . . . . . . . . . . 29
4.4.2. Recursive Puzzling Step Algorithm . . . . . . . . . . . . . . . . . 31
4.4.2.1. Example . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
4.4.2.2. Complexity . . . . . . . . . . . . . . . . . . . . . . . . . 34v

vi

Contents

4.4.3. Berry’s MRCA-Based Puzzling Step Algorithm . . . . . . . . . . .
4.5. Runtime Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.5.1. Datasets . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.5.2. Benchmark Setup . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.5.3. Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.6. Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

36
38
38
39
40
40

5. Parallelized Quartet Puzzling
43
5.1. Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
5.2. Parameter Estimation for Evolutionary Models . . . . . . . . . . . . . . . 44
5.2.1. Algorithm: Estimating Model Parameters . . . . . . . . . . . . . 44
5.3. Datasets . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
5.4. Runtime Analysis of the Sequential TREE-PUZZLE . . . . . . . . . . . . 46
5.4.1. Runtime of the Parameter Estimation . . . . . . . . . . . . . . . . 46
5.4.2. Runtime of the Quartet Puzzling Algorithm . . . . . . . . . . . . 46
5.5. Parallelizing TREE-PUZZLE . . . . . . . . . . . . . . . . . . . . . . . . 47
5.5.1. Parallelizing the Parameter Estimation . . . . . . . . . . . . . . . 47
5.5.2. Parallelizing Quartet Puzzling . . . . . . . . . . . . . . . . . . . . 47
5.6. Scheduling Algorithms . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
5.7. Efficiency of Parallel TREE-PUZZLE . . . . . . . . . . . . . . . . . . . . 50
5.7.1. Benchmark Datasets and Setup . . . . . . . . . . . . . . . . . . . 50
5.7.2. Results for Parameter Estimation . . . . . . . . . . . . . . . . . . 51
5.7.3. ML Step and Puzzling Step . . . . . . . . . . . . . . . . . . . . . 51
5.7.4. Overall Scaling of Parallel TREE-PUZZLE . . . . . . . . . . . . . 51
5.8. Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56
6. Large ML Trees from Sequences Using Quartets
57
6.1. Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
6.2. The Modified Quartet Puzzling Algorithm . . . . . . . . . . . . . . . . . 57
6.2.1. The Algorithm ModPUZZLE . . . . . . . . . . . . . . . . . . . 58
6.3. Computational Aspects . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
6.4. Some Practical Measures on the Method . . . . . . . . . . . . . . . . . . 60
6.5. Discussion and Possible Extensions . . . . . . . . . . . . . . . . . . . . . 61
7. Phylogenetic Trees from Multiple Genesets with Missing Data
63
7.1. Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63
7.2. Methods to Combine Datasets . . . . . . . . . . . . . . . . . . . . . . . . 64
7.2.1. Low Level Combination: Total Evidence . . . . . . . . . . . . . . 64
7.2.2. High Level Methods . . . . . . . . . . . . . . . . . . . . . . . . . 64
7.2.2.1. Combining Equal-Sized Datasets: Consensus Methods . 64
7.2.2.2. Combining Overlapping Datasets: Supertree Methods . . 66
7.3. Medium Level Combined Phylogenetic Analysis . . . . . . . . . . . . . . 68
7.3.1. Notation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68
7.3.2. The Combined Quartet Method to Combine Genesets . . . . . . . 70

Contents

7.3.2.1. Combining the ML Quartets . . . . . . . . . . . . . . . .
7.3.2.2. Computing the Overall Tree . . . . . . . . . . . . . . . .
7.3.2.3. Assessing Whether Genesets Can Be Combined . . . . .
7.3.2.4. Overlap-Guided Puzzling Step . . . . . . . . . . . . . . .
7.3.2.5. Relative Majority Consensus . . . . . . . . . . . . . . .
7.4. The Phylogeny of the Grasses . . . . . . . . . . . . . . . . . . . . . . . .
7.4.1. The Dataset . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7.4.2. Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7.4.2.1. Low Level Combination . . . . . . . . . . . . . . . . . .
7.4.2.2. High Level Combination . . . . . . . . . . . . . . . . . .
7.4.2.3. Medium Level Combination . . . . . . . . . . . . . . . .
7.4.3. Parameter Estimates from Poaceae Dataset . . . . . . . . . . . .
7.4.4. Combinability of the Poaceae Dataset . . . . . . . . . . . . . . . .
7.4.5. Reconstructed Poaceae Phylogeny . . . . . . . . . . . . . . . . . .
7.4.5.1. Total Evidence Tree . . . . . . . . . . . . . . . . . . . .
7.4.5.2. MRP Supertrees . . . . . . . . . . . . . . . . . . . . . .
7.4.5.3. MinCut Supertrees . . . . . . . . . . . . . . . . . . . .
7.4.5.4. Combined Quartet Tree . . . . . . . . . . . . . . . . . .
7.4.6. Other Published Poaceae Phylogenies . . . . . . . . . . . . . . . .
7.5. Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7.5.1. Problems of Dataset-Combining Methods . . . . . . . . . . . . . .
7.5.2. Comparison of the Tree Topologies . . . . . . . . . . . . . . . . .
7.5.3. Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

vii

70
71
71
72
73
73
73
74
74
74
74
76
77
78
78
78
82
82
86
88
88
89
90

8. Summary

93

A. Variables and Functions

95

B. Abbreviations

98

Bibliography101

Acknowledgments

111

List of Figuresviii

1.1. Nucleotide Database Growth . . . . . . . . . . . . . . . . . . . . . . . . .
1.2. Protein Database Growth . . . . . . . . . . . . . . . . . . . . . . . . . .

2
3

2.1.
2.2.
2.3.
2.4.

Sequence Evolution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
DNA Substitution Rates . . . . . . . . . . . . . . . . . . . . . . . . . . .
Tree Notation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Balanced and Linearized Tree Topologies . . . . . . . . . . . . . . . . . .

6
8
11
13

3.1.
3.2.
3.3.
3.4.
3.5.
3.6.

Binary Quartet Tree Topologies and Their Weights . . . . . . . . . . . .
Tree and an Induced Quartet Tree Topology . . . . . . . . . . . . . . . .
Partly and Unresolved Quartet Topologies and Their Weights . . . . . .
Likelihood Mapping . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Puzzling Step: Finding the Insertion Branch in a 4-Tree . . . . . . . . . .
Puzzling Step: Finding the Insertion Branch in a 5-Tree . . . . . . . . . .

20
20
21
22
23
23

4.1. Complexity Notations . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.2. Penalty Sums Over the Tree . . . . . . . . . . . . . . . . . . . . . . . . .

26
32

5.1.
5.2.
5.3.
5.4.
5.5.

Parallelized Workflow of TREE-PUZZLE . . . . . . . . . . . . . . . . . .
Speedup and Runtime of Parameter Estimation . . . . . . . . . . . . . .
Speedup and Runtime of the ML Step . . . . . . . . . . . . . . . . . . .
Speedup and Runtime of the Puzzling Step . . . . . . . . . . . . . . . . .
Speedup and Runtime of the TREE-PUZZLE Program . . . . . . . . . .

48
52
53
54
55

6.1. Modified Quartet Puzzling Algorithm . . . . . . . . . . . . . . . . . . . .

59

7.1. Levels of Dataset Combination . . . . . . . . . . . . . . . . . . . . . . . .
7.2. Consensus Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7.3. Adams Consensus . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7.4. MRP Coding Schemes . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7.5. MinCut Supertree Algorithm . . . . . . . . . . . . . . . . . . . . . .
7.6. ModMinCut Supertree Algorithm . . . . . . . . . . . . . . . . . . .
7.7. Overlap Graph of Gene Datasets . . . . . . . . . . . . . . . . . . . . . .
7.8. Total Evidence Phylogeny of the Poaceae Dataset . . . . . . . . . . . . .
7.9. MRP-BR Phylogeny of the Poaceae Dataset . . . . . . . . . . . . . . . .
7.10. MRP-Pu Phylogeny of the Poaceae Dataset . . . . . . . . . . . . . . . .

65
66
66
67
69
70
77
79
80
81

List of Figures

7.11. MinCut supertree of the Poaceae Dataset . . . . . . . . . . . . . . . . .
7.12. ModMinCut supertree of the Poaceae Dataset . . . . . . . . . . . . . .
7.13. Combined Quartet Phylogeny of the Poaceae Dataset . . . . . . . . . . .
7.14. Extended Majority Consensus from All Intermediate Trees . . . . . . . .
7.15. Other Poaceae Trees . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

ix

83
84
85
87
88

List of Tablesx4.1. Quartets and Topologies for an Example . . . . . . . . . . . . . . . . . .
4.2. Runtime Analysis of Puzzling Step Algorithms . . . . . . . . . . . . . . .

32
39

5.1. Runtime Profile for ML Parameter Estimation . . . . . . . . . . . . . . .
5.2. Runtime Profile for the TREE-PUZZLE Parts . . . . . . . . . . . . . . .

46
47

6.1. ModPUZZLE: Simulations with Increasing Overlap . . . . . . . . . . .

61

7.1. Sequences in the Poaceae Dataset . . . . . . . . . . . . . . . . . . . . . .
7.2. Parameters of the Poaceae Datasets . . . . . . . . . . . . . . . . . . . . .

75
76

List of Listings

4.1.
4.2.
4.3.
4.4.
4.5.
4.6.

Original Puzzling Step Algorithm . . . . . . . . . . . . . . . . . . . . . .
Split-based Puzzling Step Algorithm . . . . . . . . . . . . . . . . . . . .
Recursive Puzzling Step Algorithm . . . . . . . . . . . . . . . . . . . . .
SubtreePenalty for the Recursive Algorithm . . . . . . . . . . . . . .
MRCA-Based Puzzling Step Algorithm . . . . . . . . . . . . . . . . . . .
PropagatePenalties for the MRCA-based Algorithm . . . . . . . . .

28
30
34
35
37
38

5.1. Estimating Model Parameters in TREE-PUZZLE . . . . . . . . . . . . .

45

7.1. MinCut Supertree . . . . . . . . . . . . . . . . . . . . . . . . . . . .

69

xi

xii

List of Listings

1. Overview
1.1. Motivation
It is commonly accepted that contemporary genes, genomes, and organisms evolved fromancestors under the influence of natural selection. Consequently, the knowledge of theevolutionary tree behind their origin is crucial for understanding these entities. Evolutionary phylogenetics dealing with the reconstruction of those evolutionary relationshipsof organisms, genes, or gene families has become a basic application in many fields ofresearch.
Knowledge about the relationships within gene families plays an important role in understanding, for example, the origins of biochemical pathways, regulatory mechanisms incells as well as the development of complex systems. For example, knowing relationships
between viruses is central for understanding their ways of infection and pathogenicity.
Phylogenetic analysis is usually performed on sets of aligned related sequences fromthe primary databases like EMBL/GenBank (Stoesser et al., 2002; Benson et al., 2003,
nucleotide sequences) and SWISS-PROT/TrEMBL (Boeckmann et al., 2003, protein sequences) or databases specialized on gene families like HOVERGEN (Duret et al., 1994),
SYSTERS (Krause et al., 2002), the Ribosomal Database Project (RDP II, Maidak et al.,2001), and the European Small Subunit rRNA Database (Van de Peer et al., 2000b).
The development of efficient automated sequencing techniques and genome projects
has produced a tremendously growing amount of data in these public databases (for alist see Baxevanis, 2003) during the recent decades. While a hand-curated high-quality
database like SWISS-PROT grows about linearly (Fig. 1.2, dashed curve), other primarydatabases without the bottleneck of manual interaction like EMBL and TrEMBL show
exponential growth rates (Figs. 1.1 and 1.2).
The enormous content of the public databases serves the need for larger datasets,
since the increase of information is assumed to improve the accuracy of the results.
The datasets compiled from the different data sources can increase in two dimensions,vertically and horizontally. Vertically means that the number of taxa or species isincreased, horizontal growth describes the inclusion of more sites by obtaining either
longer sequences or by adding different sequences to each taxon. The availability oflarge datasets motivates large scale projects like the assembly of the tree of life, a global
phylogeny of all known organisms.
On the other hand, large datasets have led to severe difficulties, regarding runtime
as well as memory requirements, for the efficient phylogenetic analysis of various types

1

2

1. Overview

growth of the EMBL nucleotide database
2.5e+07
EMBL

No. Sequences

2e+07

1.5e+07

1e+07

5e+06

01/03

01/02

01/01

01/00

01/99

01/98

01/97

01/96

01/95

01/94

01/93

01/92

01/91

01/90

01/89

01/88

01/87

01/86

01/85

01/84

01/83

01/82

0

Figure 1.1.: Public nucleotide databases like the EMBL database show exponential
growth, doubling their content every nine monthof large datasets. To be able to elucidate such humongous amounts of data, it is ofutmost importance to design fast and efficient algorithms for molecular sequence analysis. This especially applies to maximum likelihood methods in phylogenetic analysis
due to their enormous need of computational resources. To make maximum likelihood
tree reconstruction available for different kinds of large datasets, this thesis introduces
four different approaches. First, the speedup of phylogenetic analysis using efficient algorithms is studied. Second, the waiting time of phylogeny reconstruction is reduced byintroducing parallel computing. Furthermore, a method is suggested that can be used
to find the coarse phylogenetic structure for a large dataset. Finally, a novel method is
proposed to combine overlapping datasets of different genes into one overall tree. All
these approaches are based on Quartet Puzzling (Strimmer and von Haeseler, 1996), aquartet-based phylogeny reconstruction method.

1.2. Organization of this Thesis
After a general introduction to biological data and phylogenetic analysis (Chapter 2)
the underlying quartet-based algorithms, namely, Likelihood Mapping (Strimmer andvon Haeseler, 1997) and Quartet Puzzling (Strimmer and von Haeseler, 1996; Strimmer
et al., 1997) are delineated (Chapter 3). This thesis approaches the phylogenetic analysisof large datasets from four different angles described in Chapters 4 to 7.
Chapter 4: Chapter 4 introduces two new algorithms which reduce the level of complexity in the puzzling step from O(n5 ) to O(n4 ). The improved algorithms are1.2. Organization of this Thesis

3

growth of public protein databases
800000TrEMBLSWISS-PROT
700000

number of sequences

600000500000400000300000200000100000

01/03

01/02

01/01

01/00

01/99

01/98

01/97

01/96

01/95

01/94

01/93

01/92

01/91

01/90

01/89

01/88

01/87

01/86

0

Figure 1.2.: While hand-curated protein databases like SWISS-PROT show only linear
growth rates, the number of protein sequences derived from the nucleotide
databases still grow exponentially (here TrEMBL, established in 1996).
compared to the original implementation and to another O(n4 ) algorithm suggested
by Vincent Berry (Strimmer, personal communication).
Chapter 5: Due to their frequently still polynomial complexity, heuristics often do not
sufficiently reduce the runtime of large dataset analysis. Hence, parallel computation has proven to be a valuable tool to decrease waiting time for the analysis.
Chapter 5 describes the parallelization of the tree reconstruction in the TREEPUZZLE package.
Chapter 6: The maximum likelihood step of Quartet Puzzling causes a major part of
the runtime of tree reconstruction. Therefore, in Chapter 6 a modified Quartet
Puzzling algorithm is suggested to decrease computation time by reducing thenumber of quartets used for phylogeny reconstruction.
Chapter 7: As mentioned, large datasets do not only consist of sequences from many
species (vertically large). Since long sequences are considered to increase information in the analysis, sequences of different genes are analyzed simultaneously(horizontally large). Yet available sequences are not spread evenly among species
of interest. Therefore, chapter 7 introduces a quartet-based method to combine
incomplete genesets for phylogenetic reconstruction. This method is discussed in
the context of total evidence and SuperTree methods.
At the end of the thesis the results will be summarized and a short outlook to future
work will be given (Chapter 8).

4

1. Overview

2. Molecular Phylogenetics: A General
Introduction
This chapter will present some basics of phylogenetic analysis. In the first section biological sequence data, the concept of molecular evolution, and a way to model theevolutionary process will be introduced. The basic notations on trees will be defined in
the second part of the chapter, followed by the explanation of different types of phylogenetic methods, and a brief introduction to maximum likelihood phylogenies. Finally,some facts about the complexity of phylogenetic reconstruction will be given.

2.1. Biological Sequences and Molecular Evolution
2.1.1. Types of Biological Data
To study evolutionary relationships, very diverse biological data are used like morphological characters, binary data, genomic gene order, nucleotide, and protein sequencedata.
The first data used were morphological characters of the species under interest. Morphological characters have the advantage to be easily obtainable by eye or by microscopy
without molecular laboratory work.
In the early times of molecular phylogenetics, coarse grain genetic data like binary
characters of presence and absence of restriction sites played an important role. Veryrecently, gene order data has been used to reconstruct phylogenetic relationships fromthe order of genome rearrangements, so-called breakpoint phylogenies (Blanchette et al.,
1997; Sankoff and Blanchette, 1998).
Today the vast majority of data in phylogenetic reconstruction is biological sequence
data like nucleotide or protein sequences as stored in the public databases (see p. 1). One
major advantage of biological sequence data is the increase of phylogenetically relevant
information due to the large number of characters, that is, sites or residues, that can beobtained by sequencing.
Molecular sequences are coded as strings of literals from an alphabet A representing the order in which the building blocks are connected in the molecules. Nucleotide
sequence data (DNA and RNA) are coded with an alphabet of the four nucleotides
Adenine, Guanine, Cytosine, and Thymine in DNA or Uracil in RNA. Each nucleotide

52. Molecular Phylogenetics: A General Introductiont... A G C C

T ...

s

... A G C A

T ...

1

1

t

2

s

2

... G G C A

n

su

bs

tio

le

tit
u

time

de

tio

n

t

s

0

0

substitution

6

T ...

... A G−A

T ... s’2

Figure 2.1.: An example for sequence evolution. Sequence s0 evolves to s1 in time t1 −t0 .
Then s1 evolves divergently to s2 and s02 after a speciation/duplication eventin time t2 − t1 .
or base is denoted by its first letters, A = {A, G, C, T /U } (Liébecq, 1992). Nowadays
nucleotide sequences can be efficiently obtained using fast (automated) sequencing techniques (Sanger et al., 1977 and, e.g., Ansorge et al., 1993).
Protein data is encoded over an alphabet of 20 amino acids (Liébecq, 1992). Eachamino acid is represented as a one-letter or as a three-letter literal (’A’ or ’Ala’ for
Alanine, ’P’ or ’Pro’ for Proline etc.). Although there has been a lot of work on protein
sequencing (Edman, 1950; Sanger and Thompson, 1953) and its automation (Edman and
Begg, 1967) since the 1950s, the indirect strategy of sequencing protein coding genes withthe method of Sanger et al. (1977) and the subsequent translation into amino acids has
proven to be faster and has become the usual way of deriving protein sequences.

2.1.2. Sequence Evolution
Before the division of any cell, plastid, or mitochondrion its genome has to be replicatedto be inherited to the daughter cells or organelles. In spite of a proof reading machinery
the process of copying is not error-free. Additionally, damages are introduced to the
DNA by mutagens such as certain chemicals or UV light. Hence, the genomic sequence
accumulates mutations as traces of evolutionary development. Although more complex
mutations like rearrangements, duplications, and inversions are possible on chromosome
level, only point mutations affecting single spots of the DNA are commonly considered.
Point mutations aresubstitutions – the exchange of one character state by another (Fig. 2.1)
insertions– inserting one or more charactersdeletions– deleting one or more character usually indicated by gap characters ’-’in the sequence. (Fig. 2.1)

2.1. Biological Sequences and Molecular Evolution

7

Insertions and deletions cannot be distinguished when examining two contemporary
sequences (e.g., s2 and s02 in Fig. 2.1). To decide whether an insertion or deletion
has happened, one needs information on the sequence of the common ancestor (s1 in
Fig. 2.1) which, however, is usually unknown. For that reason insertions and del etion are
generalized as indels. Indels and substitutions generally apply to nucleotide as well asprotein data. Although mutations happen at the genomic level, they can affect the aminoacid sequence translated from protein coding genes. Nucleotide substitutions in protein
coding genes that do not result in an amino acid substitution are called synonymous orsilent mutations in contrast to non-synonymous mutations.

2.1.3. Multiple Sequence Alignment
Homologous sequences, that means, sequences related by a common ancestor sequence,
are typically presented in a multiple sequence alignment (MSA). An alignment is a data
matrix in which homologous characters, also called sites, of the sequences are aligned in
the same column. For instance, the alignment of the sequences from Fig. 2.1 is
s0 :
s1 :
s2 :
s02 :

...
...
...
...

A
A
G
A

G
G
G
G

C
C
C
−

C
A
A
A

T
T
T
T

...
...
...
...

(2.1)

To present indels, an additional gap character ’-’ is used. Sequence alignments serve
as input to almost every sequence analysis and programs like CLUSTALW (Thompson
et al., 1994) or DIALIGN (Morgenstern, 1999) are available to compute an MSA for acollection of sequences. However, there is still no standardized approach how to treatindels in the subsequent analysis. Typically columns with gaps are either discarded
from the alignment, or indels are considered as ’wildcard’ characters substituted by a
distribution of character states. In the following we will for simplicity delete sites withindels.

2.1.4. Modeling Molecular Evolution
To reconstruct evolution, the evolutionary process is modeled as an evolutionary Markov
process (EMP, Tavare, 1986; Müller and Vingron, 2000), also known as substitution
model. An EMP is a time homogeneous, calibrated, stationary, reversible Markov process.
It is assumed that the frequencies πρ of the characters ρ ∈ A are in equilibrium
and remained stationary during evolution and among the sequences (stationarity). The
vector π = (π1 , π2 , . . . , π|A| ) of the character frequencies is called stationary distribution.

8

2. Molecular Phylogenetics: A General Introduction
H
N

N
P

N

H

P

A

N

N
S

γ

N

S

O

N

G

N

H

N

H

N

3’

H

3’

β
N

H

CH 3

C
N
P

σ

α

τ

H

ρ

T
N
P

N
S

O

O

S

3’

H

N
O

3’

Figure 2.2.: DNA substitution rates
Since we infer evolution backward in time and generally have no means to deriveancestral sequences, mutations are modeled to be reversible, that is, mutations from ρto σ are assumed to be equally likely as the back-mutation from σ to ρ (reversibility)πρ Pρσ (t) = Pσρ (t)πσρ, σ ∈ A.

(2.2)

This assumption is also known as detailed balance.
Since evolution is modeled as a Markov process, it evolves without memory, that is,for the development from sequence s1 to s2 in Fig. 2.1 the character state of s0 does notmatter, only the state s1 at t1 is of importance. Furthermore, it is assumed that theMarkov process acts on each character of the sequence independently from the others(independence).
The following applies to DNA as well as protein evolution but for simplicity we consider
the smaller nucleotide alphabet A = {A, C, G, T } only. The probability of substitutionof any character ρ by character σ in evolutionary time t is described by a transition
probability matrix P . Given the stationary distribution π and substitution rate matrix
R (see Fig. 2.2 for substitution rates)

 A
−
 β
R= 
 γ
α

C
β
−
τ
ρ

G
γ
τ
−
σ

T 
A
α

ρ  C
σ  G
−
T

(2.3)

2.1. Biological Sequences and Molecular Evolution

9

an instantaneous rate matrix Q can be constructed such that


Q=


−

P

k6=A QAkµRCA πAµRGA πAµRT A πA

µRAC πCµRAG πGµRAT πTPµRCG πGµRCT πT − k6=C QCkPµRGC πC− k6=T QGkµRGT πT PµRT C πCµRT G πG
− k6=T QT k

(2.4)

Hence, for the row sums the following applies.
X

for each ρ ∈ A.

Qρσ = 0

(2.5)

σ∈A

The factor µ in Eq. 2.4 is calibrated such thatXQρσ = −XQρσ = 1ρ, σ ∈ A.

(2.6)

ρ

ρ6=σ

The transition probability matrix P is constructed from Q
P (t) = eQt .

(2.7)

Sometimes (cf. Müller and Vingron, 2000) µ is calibrated such thatXπρ Pρσ (1) = 0.01

ρ, σ ∈ A.

(2.8)

ρ6=σ

That means, evolutionary time is calibrated such that in time t = 1 (1 PAM) a sequenceof length m accumulates m/100 (= 1%) mutations.
P also has the following properties:
Xπρ Pρσ (t) = 1.0

for each t ≤ 0 and ρ ∈ A

(2.9)

σ∈A

and

lim Pρσ (t) =t&0

0.0, if ρ 6= σ1.0, if ρ = σ

ρ, σ ∈ A.

(2.10)

Note, that an EMP is uniquely described by Q or P (t), since π follows the equations
πQ = 0 and πP (t) = π for all t > 0 (Tavare, 1986).
As mentioned above, Eq. 2.2 to Eq. 2.10 also apply for amino acids, but 20 × 20matrices are used instead of 4 × 4-matrices, due to the larger alphabet size |A|.
A number of models has been suggested both for nucleotide as well as amino acidevolution. The models vary in complexity regarding their parameters. Common DNA

10

2. Molecular Phylogenetics: A General Introduction

models (listed in ascending complexity) are JC69 (Jukes and Cantor, 1969), Kimura-2Parameter model (K2P, Kimura, 1980), HKY85 (Hasegawa et al., 1985), TN93 (Tamura
and Nei, 1993), and the General Time Reversible model (GTR, Rodriguez et al., 1990)
as used in Eq. 2.3 and Fig. 2.2.
To model amino acid evolution, matrices are designed for different scenarios. While
the Dayhoff model (Dayhoff et al., 1978) and JTT (Jones et al., 1992) are generally applicable, VT (Müller and Vingron, 2000) and WAG (Whelan and Goldman, 2001) model
distant evolution. Specialized matrices like mtREV (Adachi and Hasegawa, 1996a) and
cpREV (Adachi et al., 2000) model mitochondrial and chloroplast evolution respectively.
Transition probability matrices are applied to database searches and alignment construction. Furthermore, they are used to compute evolutionary distances between sequences and to evaluate trees for phylogenetic reconstruction in a maximum likelihood(ML) framework. For more detailed information on transition probability matrices see,e.g., Ewens and Grant (2001), Graur and Li (2000), Durbin et al. (1998), and Tavare(1986).
Additional to the transition probability matrix, other assumptions can be made tomore realistically model the evolutionary process. It is known that sequence positions
may evolve with different rates. To account for this fact, each site in an alignment may
also get its own rate specific factor. Typically one assumes that rates are distributed
according to a Γ-distribution, where one parameter α describes the amount of rate
heterogeneity (Gu et al., 1995; Yang and Kumar, 1994). Apart from that, methods
have been suggested to estimate site-specific rates for each alignment column without
assuming a Γ-distribution (e.g., Van de Peer et al., 2000a; Meyer and von Haeseler, 2003,
and references therein).

2.2. Phylogenetic Tree Reconstruction
Recent taxa and their genetic sequences evolved from ancestors. Hence, tree diagrams– phylogenetic trees – are a logical way to present evolutionary relationships. Before
going into detail how phylogenetic trees are reconstructed, the basic notations on trees
will be defined.

2.2.1. Notation on Phylogenetic Trees
The basic tree notations used throughout this thesis will be defined by means of theexamples in Fig. 2.3. Certain characteristics of tree topologies will be explained using
the trees in Fig. 2.4. If special notations are needed later on they will be defined in thecorresponding chapters.

2.2. Phylogenetic Tree Reconstructiona11bleaf/taxon

A

B

C

D

E

F

G

leaves/taxa

G

F

E

D

A

B

C

4

external
branch1bifurcation

2

4

internal branch

3

2

3

inner node

inner node

multifurcationrootcinternal
branchC

root

multifurcationdDexternal
branch inner
node

H
E

1

B

2

bifurcation

3
leaves/taxa

external
branch1

inner node4

leaves/taxa

5

F

A

I
external
branchG

J

6K
internal
branchFigure 2.3.: Tree notation on different rooted tree diagrams (a,b), an unrooted tree (c),and unrooted quartet tree (d)Trees, Nodes, and Edges A tree is an acyclic graph T = (V, E) consisting of a set ofnodes V and a set E of edges linking two nodes each. If a trees has a root node ∈ V it iscalled a rooted tree represented by a directed acyclic graph (DAG) for which all edgese ∈ E are directed away from the root. Edges are denoted by their adjacent nodes, e.g.,e(1,2) describes the edge linking nodes 1 and 2. In a rooted tree e(•,2) describes the edge
ending at node 2.
The degree of a node v is defined as the number of edges connected to v. External
nodes, also called leaves, are connected to one branch only, they have degree 1. Nodeswith degree > 1 are internal nodes. Except for the root node in rooted trees all inner
nodes will have at least degree 3. A node with three adjacent branches (degree 3) iscalled a bifurcation. Nodes with degrees > 3 are called multifurcations (Fig. 2.3). Atree where all inner nodes have degree 3 is called a bifurcating or binary tree. Otherwise
it is a multifurcating tree which is not completely resolved.
Evolutionary Trees, Leafsets, and Leaf-Labeled Trees In evolutionary trees, also
called phylogenies, the nodes correspond to taxa and their ancestors. The edges, also
called branches, depict the development between two nodes. Edges can be assigned
weights to express evolutionary time. Evolutionary time is displayed as branch lengths(Fig. 2.3a,c). In Fig. 2.3b, only the vertical branch lengths depict time, while thehorizontal part does not count.
Generally, sequences represent the taxa they are derived from. The leaves of the12

2. Molecular Phylogenetics: A General Introduction

phylogenetic tree are labeled by contemporary sequences (A to G in Fig. 2.3), while theinternal nodes (nodes 1 to 4) correspond to the hypothetical ancestor sequences. Since
we have no knowledge about ancestral sequences, evolutionary trees are leaf-labeledtrees. Hence, the terms taxon, species, sequence, leaf, and external node are usedinterchangeably.
To access leafsets of trees we define the following. T (S) is a tree with leafset S, andL(T ) describes the leafset of tree T . Hence, L (T (S)) is exactly S.
A tree with n leaves is called an n-tree. The topology of an n-tree is denoted by Tn .
Accordingly, the size of the leafset |L(Tn )| of an n-tree is n.
An n-tree has a maximal number of n − 2 internal nodes and 2n − 3 edges if it iscompletely resolved. A completely unresolved tree comprises the minimal number of
one internal node and n edges, i.e., it has only external edges.

Biological Root and Outgroup Rooting Since the sequence at the general ancestor is
unknown the placement of the root often remains a problem. Therefore, the phylogeny
reconstruction programs usually produce ’unrooted trees’ like in Fig. 2.3c. If one wants
to determine the order of evolutionary events additional information is necessary to place
the root on the tree. One approach is outgroup rooting. Outgroups comprise one or moretaxa that are known to have branched off the tree of life before the interesting group of
taxa in the analysis. The root is then placed on the branch between the outgroup andthe rest of the tree. The unrooted tree in Fig. 2.3c has been rooted by the outgroupsequences A, B, and C in the trees in Fig. 2.3a and b.

Splits and Bipartitions The internal edges induce bipartitions on the leafset, thatmeans, when cut they split the set of leaves into to two disjoint sets. In Fig. 2.3, theedge e(1,2) between nodes 1 and 2 splits the leaf sets of the trees into the disjoint subsets{B, C} and {A, D, E, F, G}. This bipartition or split is denoted as BC|ADEF G. In
case of trees of four sequences this notation exactly describes the tree topology since
there is only one internal branch possible (e.g., HI|JK in Fig. 2.3d). To describe thesubsets induced by a splitting edge e, we define
e
e
r e
r e

=
=
=
=

{taxa, on the left of edge e}{taxa, on the right of edge e}{taxa, on the same side as the root of edge e}{taxa, not on the same side as the root of edge e}(2.11)(2.12)(2.13)(2.14)

Although the sets e and e are not uniquely defined, the two notations are sufficient todescribe the two set of leaves induced by an edge e.

2.2. Phylogenetic Tree Reconstruction
a

13
b

A

B

C

D

E

F

G

H

c

A

B

C

D

E

F

G

H

d
A

E

B

F

C

G
D

A

H

B

G
C

D

E

F

H

Figure 2.4.: Balanced and linearized tree topologies: balanced rooted tree (a), linearized
rooted tree (b), balanced unrooted tree (c), linearized unrooted tree (d) also
called caterpillar or comb tree
Rooted Trees, Ancestors, and Clusters In a rooted tree, succ(w) denotes the set of
immediate descendants of node w, e.g., succ(3) = {4, E, F } in Fig. 2.3a and b.
Nodes in rooted trees allow the notation of common ancestors. The most recent
common ancestor (MRCA) of a set of leaves or nodes is defined as the ancestral node
furthest from the root that has all these leaves as descendants. In Fig. 2.3, the mostrecent common ancestor of the nodes D, 4, F is node 3, i.e., MRCA(D, 4, F ) = 3, while
MRCA(B, E) = root.
The set of all descendants of an ancestral node is called a monophyletic group, orclade. Taxa A, B, and C are forming a monophyletic group while D, F , and G are notmonophyletic. The set of descendant leaves of an internal node is also called the clusterof that node. The cluster of node 1 in Fig. 2.3 is C n (1) = {B, C}. Since in rooted trees
edges are directional, we also speak about the cluster of an edge, e.g., the cluster of theedge e(2,1) is also C e (e(2,1) ) = {B, C} which is identical to r e(2,1) .
Binary Tree Topologies, Diameter, and Height Two topologies among the binary
trees are extreme with respect to the diameter, the diameter of a tree being defined as
the longest path between any two nodes in the tree. First, there are fully balanced trees(Fig. 2.4a,c) where all leaves have the same height. The height h is defined as the numberof edges between that leaf and the root in rooted trees or that leaf and the central edge
in unrooted trees. Second, there are linearized trees also called caterpillar or comb-like14

2. Molecular Phylogenetics: A General Introduction

trees, which are highly unbalanced (Fig. 2.4b,d). All other binary tree topologies areintermediates between these extremes. The number of binary unrooted n-trees B(n) is

B(n) =nY(2k − 5) =k=3

(2n − 5)!
2n−3 (n − 3)!

(2.15)(Felsenstein, 1978). The number of possible rooted n-trees is B(n + 1), considering theroot as an additional leaf.

2.2.2. Types of Phylogenetic Methods
A rich variety of tree reconstruction methods based on sequences has been developed(for an overview see Swofford et al., 1996; Baxevanis et al., 2002), which fall into three
categories, (a) maximum parsimony methods (MP), (b) evolutionary distances (ED)between sequences, and (c) approaches applying the maximum likelihood (ML) principle.
MP tries to find the tree which explains the data with the least mutations. ED methods
calculate pairwise distances between the taxa and construct a tree from these distances.
ML aims to find the tree that gains the maximum likelihood to have produced theunderlying data. Most of these methods use objective functions to evaluate trees andaim to find trees optimal in respect to the objective functions.
In the following we will describe how trees a evaluated using maximum likelihood. ML
analysis is statistically well founded and used in many fields of research (e.g., Cramer,
1989; Eliason, 1993; Mendel and Burrus, 1990; Ewens and Grant, 2001). Since ML has
been made applicable for reconstructing phylogenies of molecular sequences in 1981 byFelsenstein, it has shown to give accurate results in practice and to be quite robust
against several sources of error. Although phylogenetic methods applying ML are computationally very expensive, their use in phylogeny reconstruction has increased. They
are a routine method in phylogenetic analysis.

2.2.3. Maximum Likelihood on Phylogenetic Trees
In this section ML analysis is described briefly, for more details on maximum likelihood
phylogenies refer to Felsenstein (1981), Goldman (1990), and Swofford et al. (1996).
In ML analysis a likelihood value is used to describe the goodness of fit of the hypothesis H and the data A. The likelihood of the hypothesis H given the data A is definedas the probability of the data A which in turn is the product of the partial probabilities
for each data sample Ac (Goldman, 1990)
L(H|A) = Pr(A|H) =Yc

Pr(Ac |H).

(2.16)

2.2. Phylogenetic Tree Reconstruction

15

The basic principle of ML estimation is that the best estimate of an unknown parameter of the hypothesis H is the one making the analyzed data most probable, whichmeans, its likelihood value is maximized.
As data in molecular phylogenetics serves an alignment A comprising m columns of aset S = {s1 , s2 , . . . sn } of n aligned homologous sequences like

s1
s2
..
.
sn

A1 . . . Ac−1 Ac Ac+1
A ... G
T
A
A ... G
C
T
..
..
..
..
.
.
.
.
T ... GCG

. . . Am
... C
... C
..
.
...

(2.17)C

To compute the likelihood, we need a hypothesis H of the process how the data A
came into existence. We assume that the contemporary sequences evolved according toa phylogenetic tree T under the influence of an EMP M (2.1.4). The phylogenetic tree T
consists of an n-tree topology and a set of branch lengths t. tvw denotes the length of the
branch e(v,w) that connects the two nodes v and w. Each of the m alignment columns Ac
is regarded to be an independent sample from the EMP acting along the tree. Although
the hypothesis H includes both, the tree T and the EMP M , for simplicity we willconsider the tree T only, we assume M to be fixed.
To evaluate T , its likelihood value is calculated. If the tree T (topology, branch
lengths) and the character states at all nodes are known, the calculation is straightforward. For each data sample the likelihood is computed by multiplying the probabilitiesPρσ (tv,w ) (Eq. 2.7) for all edges e(v,w) , where ρ and σ are the character states at theadjacent nodes v and w and tv,w is the evolutionary time between these node. Finally,
the likelihoods of all data samples are multiplied to gain the overall likelihood of thetree.
However, since we do not know the ancestral sequences we have to sum over all
possible character states at the internal nodes. The possible states are the elements of
the alphabet A, that is, 20 amino acids for protein data or 4 nucleotides for DNA.
The likelihood LT = L(T |A) of the tree T is computed from the alignment A using
the probability matrix P (te ) for each edge e (Eq. 2.7) and the stationary distribution
π of M as described by Felsenstein (1981). LT is computed recursively in a post-order
traversal of the tree, i.e., the probabilities of the states at the nodes will be computed
after calculating the probabilities at the descendant nodes. The recursive calculation
will proceed as follows.
Assume T a rooted n-tree. At the moment we only consider one column Ac of thealignment, which is possible under the assumption that the sites are independent. Ac,kdescribes the character state of sequence sk in that column.

16

2. Molecular Phylogenetics: A General Introduction

Assuming node v to be a leaf representing sequence sk , the likelihood Lv (ρ) of thecharacter state ρ at leaf v is determined from the alignment column since its state is
known.

1 if ρ = Ac,kLv (ρ) =ρ ∈ A.
(2.18)
0 otherwise
If node v is an inner node, the probability of state ρ at node v having given rise to statesσ ∈ A at a descending node u is computed by summing up the products of probabilityPρσ (tv,u ) that ρ is σ after time tv,u and the partial likelihood Lu (σ) to find state σ atnode uX
Lv (ρ) =Pρσ (tv,u )Lu (σ).
(2.19)σ∈A

Since an inner node v has at least two descendants, the likelihood of any state ρ atnode v is calculated by multiplying all probabilities for all succeeding nodes w ∈ succ(v)!
YX
Lv (ρ) =Pρσ (tv,w )Lw (σ) .
(2.20)w∈succ(v)σ∈Ac
The likelihood LA
T = Pr(Ac |T, M ) of the alignment column Ac to be produced along
tree T by an EMP M is the sum of the partial likelihoods of the states σ ∈ A
X
c
LA
=
Pr(A
|T,
M
)
=
πσ Lroot (σ).
(2.21)cTσ∈A

From the likelihoods of the alignment columns, one computes the overall likelihood valueof A for the tree T
LT = Pr(A|T, M ) =

m
Y

k
LA
T

(2.22)

k=1

or the log-likelihood`T = ln LT =

m
X

k
ln LA
T .

(2.23)

k=1

This objective function enables the optimization of the parameters of the hypothesis.
The branch lengths are optimized to find the maximum likelihood value of tree T under
the assumption of a fixed model of evolution. Since optimizing all branch lengths at once
is not feasible, Felsenstein (1981) suggested to optimize the branches one by one. This
iterative task is accomplished by numerical methods such as Expectation Maximization(EM, Dempster et al., 1977; Felsenstein, 1981) or Newton-Raphson-optimization (Press
et al., 1988; Olsen et al., 1994) applying Felsenstein’s (1981) Pulley Principle.

2.2. Phylogenetic Tree Reconstruction

17

2.2.4. Complexity of Phylogenetic Analysis
Using an objective function like the one described in the previous section, a given tree
can easily be evaluated. However, to find the optimal tree gaining the highest likelihood,
one would have to examine all possible n-trees. As described in Eq. 2.15, the numberof bifurcating unrooted trees B(n) can be very large. For 10 sequences more than 2
million trees exist. For 55 taxa the trees outnumber the estimate of 1081 atoms in the
known universe. Considering multifurcations makes the growth rate even more dramatic.
Obviously, exhaustive tree searches are applicable for small data sets only.
To make things worse, the problem of finding an optimal phylogeny has been shown
to be NP-complete for quite a number of reconstruction methods (Graham and Foulds,
1982; Foulds and Graham, 1982; Day and Sankoff, 1986; Day, 1987). NP-complete
problems are those for which an efficient, i.e., polynomial-time algorithm is assumed to
be impossible. In other words, it might be necessary to examine exponentially many
trees to find the optimum. This should, in principle, also be true for ML methods.
In order to reduce the computational burden and to limit the vast number of trees
to be examined heuristics have been suggested to search the tree space. Beside others a number of ML-based heuristics have been introduced such as stepwise insertion
with local and global optimizations (DNAML, Felsenstein, 1981), the Quartet Puzzling
algorithm (TREE-PUZZLE, Strimmer and von Haeseler, 1996; Schmidt et al., 2002),
and star decomposition (MOLPHY; Adachi and Hasegawa, 1996b). Recently, Bayesian
approaches (MrBAYES, Huelsenbeck and Ronquist, 2001), genetic algorithms (Lewis,1998), and simulated annealing (Salter and Pearl, 2001) have entered the field. Note,that heuristics cannot guarantee to find the best tree.

18

2. Molecular Phylogenetics: A General Introduction

3. The Quartet Puzzling Algorithm
Since the novel methods and extensions that will be introduced throughout this thesis (chapters 4–7) are based on the Quartet Puzzling algorithm, the relevant methods
implemented in the TREE-PUZZLE package will be described in this chapter. Firstthe relevant notations on quartets are presented. Thereafter, the methods Likelihood
Mapping and Quartet Puzzling will be explained.

3.1. Introduction
While most tree reconstruction methods aim at reconstructing the overall tree for nsequences by optimizing a global objective function, quartet methods break down thistask into smaller parts. The methods attempt to reduce the problem of evaluatingcomplex objective functions (cf. 2.2.3) for a lot of large trees to the computation ofquartet trees, i.e. four-species trees, and to build the overall tree from the collection ofquartet trees. The quartet tree topologies serve as building blocks of larger trees.
Let S = {s1 , . . . , sn } be the set of n homologous (i.e. related) aligned sequences. Q
denotes the set containing all n4 quartets of S. For each quartet q = {a, b, c, d} ∈ Qthree binary tree topologies exist (Fig. 3.1). A quartet topology is uniquely describedby the bipartition induced by the internal edge. Hence, quartet topologies are denotedby ab|cd, ac|bd, ad|bc (Fig. 3.1), where a, b, c, d ∈ S are the labeled leaves of the trees.
Given a tree of n sequences, there are n4 different quartets. For each quartet q ={a, b, c, d} derived from an n-tree there is an induced unique quartet topology (Fig. 3.2).
nFrom the set of 4 induced quartet topologies, the full tree can be unambiguously
reconstructed (Bandelt and Dress, 1986). Since the overall tree is not known the set ofquartet topologies has to be inferred from aligned sequences. Several criteria exist tofind the set of quartet tree topologies supported by the input dataset. Possible methods
are again evolutionary distances, the parsimony criterion, or the maximum likelihood
framework.
If each constructed quartet tree agrees with the unknown induced quartet tree, thenreconstructing the n-tree is straight–forward. However, the stochastic process of evolution makes it unlikely to reconstruct each quartet correctly. This, on the one hand,leads to the reconstruction of contradicting quartet topologies. Consider, for instance,
AB|CD, AB|CE, and AE|BD to be reconstructed quartet trees from sequence data.
While the first two support a split AB|CDE, this is contradicted by the placement oftaxa B and E in the third tree. Thus, it is not possible to unambiguously join these

19

20

3. The Quartet Puzzling Algorithm

ac|bd

ab|cd

ad|bc

a

c

a

b

a

b

b

d

c

d

d

c

(1, 0, 0)(0, 1, 0)(0, 0, 1)

Figure 3.1.: The three different binary (fully resolved) tree-topologies for the quartet
q = {a, b, c, d}. The generic weights (wab|cd , wac|bd , wad|bc ) used for the leastsquare fit are given below the topologies.
A

J

E

I

B

A

I

D

G

H
C

F

G

D

Figure 3.2.: An example tree with 10 leaves, and the induced topology for a quartet{A, D, G, I}.
trees. On the other hand, it might often not be possible to unequivocally resolve which
topology is supported by the input data. Two or even all three topologies might gain
almost equal support by the data. These so-called unresolved cases can be depicted by
network or star tree topologies representing the neighborhood relations from more than
one binary tree topology (Fig. 3.3). Most quartet methods, however, take into accountonly one best topology for each quartet.
A number of quartet-based tree reconstruction methods have been developed to construct trees from sets of 4-trees (e.g., Sattath and Tversky, 1977; Fitch, 1981; Bandelt
and Dress, 1986; Dress et al., 1986; Strimmer and von Haeseler, 1996; Ben-dor et al.,
1998; Ranwez and Gascuel, 2001). The developments and enhancements presented inthis thesis are based on the quartet methods implemented in the TREE-PUZZLE package (Strimmer and von Haeseler, 1996; Schmidt et al., 2002; Schmidt and von Haeseler,
2003).

3.2. Methods from the TREE-PUZZLE Package
3.2.1. Likelihood Mapping
While distance-based methods or the parsimony criterion can be used to decide which of
the three informative tree topologies is best supported for one quartet, TREE-PUZZLE
uses the maximum likelihood framework to select.

3.2. Methods from the TREE-PUZZLE Package

ab|cd
ac|bd
a

b

ab|cd
ad|bc
c

a

d

b

( 12 , 12 , 0)

( 12 , 0, 12 )

21

ab|cd
ac|bd
ad|bc

ac|bd
ad|bc
d

a

c

c

d

b

(0, 12 , 12 )

a

c

b

d

( 13 , 13 , 13 )Figure 3.3.: The three different partially resolved tree-topologies and the completely unresolved one for the quartet q = {a, b, c, d}. The supported topologies combined in the respective scenario are given above, the generic weights (wab|cd ,wac|bd , wad|bc ) used for the least-square fit are given below the topologies.
For a quartet q = {a, b, c, d} the trees ab|cd, ac|bd, and ad|bc are evaluated. That is,
maximum likelihood values
Lab|cd , Lac|bd , Lad|bc(3.1)

or the corresponding log-likelihoods`ab|cd , `ac|bd , `ad|bc(3.2)

are computed for each topology according to 2.2.3.
To take into account that two or even all three topologies can receive almost identicalsupport, Strimmer et al. (1997) computed the Bayesian weights or posterior probabilitiesof the likelihoodspτ =LτLab|cd + Lac|bd + Lad|bcτ ∈ {ab|cd, ac|bd, ad|bc}.

(3.3)

For equally supported quartet trees they decided randomly which topologies to choose.
The rational behind this strategy is as follows. If one topology τ is clearly supported its
weight is pτ ≈ 1.0 while the other two are approximately zero. Such a quartet is called
fully resolved. In the case of a partly resolved quartet, that is, two topologies τ and τ 0
are equally well supported their weights are pτ ≈ pτ 0 ≈ 0.5. If all three topologies get
almost equal support (unresolved quartet) each of the weights will be approximately one
third.
Accordingly, for the 7 possible scenarios (3 three fully resolved, 3 partly resolved, and
1 unresolved state; cf. Figs. 3.1 and 3.3) a set of generic Bayesian weights (wab|cd , wac|bd ,wad|bc ) are used to discriminate the scenarios. The generic weights are chosen 1, 12 , or 13
as shown in Figs. 3.1 and 3.3. The supported quartet topologies are selected according
to the scenario with the least square distanceXD=(pτ − wτ )2(3.4)
τ ∈{ab|cd,ac|bd,ad|bc}

22

3. The Quartet Puzzling Algorithm

p =
1.0
1

1

L1
L 1 + L 2+ L 3

a

c

b

d

a

b

a

b

d

c

c

d

0.33
1.0

0.33

p =
2

0.33

1.0

p =
3

L2
L 1 + L 2+ L 3

a

b

d

c

a

b

d

c

a

c

d

b

a

b

c

d

L3
L 1 + L 2+ L 3

Figure 3.4.: How likelihood weights are plotted in a likelihood mapping diagram. Leftside: likelihood weight plotted in a three-dimensional coordinate system.
Right side: the simplex and its areas and the corresponding quartet topologies. The gray triangles are identical, only viewed from different angles.
between the generic weights wτ and the posterior probabilities pτ .
The posterior probabilities from Eq. 3.3 have the nice property of summing up to
1.0. This gave rise to a method called likelihood mapping (Strimmer and von Haeseler,1997) that visualizes the quartet likelihoods. If three values which add up to 1.0 are
plotted into a 3-dimensional coordinate system, all these points fall exactly into thetriangle spanning between the coordinates (1, 0, 0), (0, 1, 0), and (0, 0, 1). This triangle,
also called a simplex, can be printed directly by changing the angle (see Fig. 3.4). The
seven areas of the likelihood mapping diagram present the seven scenarios. Points inthe corners show unambiguously supported quartet topologies. The rectangular areas
display the partly resolved quartets, while the center comprises all unresolved quartets.
By plotting all or a large random subset of the quartets, likelihood mapping visualizes
the phylogenetic information of a dataset. The more points fall into the middle ofthe diagram, the less suitable is the dataset for phylogenetic analysis. Moreover, byclustering the sequences into 2 to 4 clusters, the relationships and monophyly of theseclusters can be examined. Each of the corners is assigned a certain cluster relationship,
and the number of dots in that area shows the respective quartet support.

3.2.2. Quartet Puzzling
The tree reconstruction method, implemented in TREE-PUZZLE, first reads a multiple
sequence alignment and estimates missing parameters according to the user specified
model of evolution. The parameter estimation process will be described in chapter 5.
The actual tree reconstruction, the Quartet Puzzling algorithm, comprises three steps:(1) maximum likelihood step (ML step), (2) puzzling step, and (3) consensus step (Strimmer and von Haeseler, 1996; Strimmer et al., 1997).

3.2. Methods from the TREE-PUZZLE Package
a)

A

23

b)

C

c)

A

C

AE | BC

1

1

A

C

AE | BD

1

B

D

1

2
2

B

D

1

B

D

AB | CD
AC | DE

A

E

A

C

C
1

4

3

B

2

f)

C
1

2

B

D

A
BD | CE

2

3

2

D

1

B

D

e)

d)

Figure 3.5.: Evaluation of the best insertion branch by penalizing. In this example
AE|BC, AE|BD, AC|DE, and BD|CE are the supported quartet topologies found in the previous step evaluating the quartets. The penalty neighbors are marked bold-face.
F

B
A

C

E

D

AC | BF
AD | BF
AE | BF
AF | CD
AE | CFAE | DFBC | EF
BF | CD
BF | DECD | EF

B
1
A

5

3
4

E

4

5

C

5
D

Figure 3.6.: Evaluation of the best insertion branch in a 5-tree by penalizing. The penalty
neighbors in the quartet topologies are marked bold-face.
3.2.2.1. ML Step
In the ML step, the set of quartet topologies is computed which act as building blocks
to reconstruct the overall tree topology. To do so, for each possible quartets q ∈ Q thethree binary topologies are evaluated by ML (cf. 2.2.3). The supported topologies arechosen applying equations Eq. 3.3 and Eq. 3.4 (cf. 3.2.1). As described before, thisapproach also selects the second best or even all three quartet topologies if there is noclear support for a single topology (Strimmer et al., 1997).
3.2.2.2. Puzzling Step
Starting from a supported topology of the first quartet {s1 , s2 , s3 , s4 }, sequences areadded one after another as follows (Strimmer and von Haeseler, 1996; Strimmer et al.,1997).

24

3. The Quartet Puzzling Algorithm

Assume X the next sequence to be added to the tree Ti that already contains the firsti sequences. For all quartets of the typeqX ∈ { {a, b, c, X} | a, b, c ∈ L(Ti ), X to be inserted} = Q(X,Ti )(3.5)

the supported topology is used to determine the optimal insertion point for sequence
X (see Figs. 3.5 and 3.6 for examples). Assume ab|cX the supported topology of aqX ∈ Q(X,Ti ) . Then for each such quartet qX ∈ Q(X,Ti ) a penalty of 1 is added to everyedge on the path between the two taxa a and b, which are not neighbors of X in the
supported topology. This path between a and b is called penalty path, while a and b
are called penalty neighbors. If more than one topology is supported for a quartet, one
topology is chosen randomly. The sequence X is added to an randomly chosen edgee ∈ Emin , where Emin is the set of edges with minimum penalty. The procedure ofsequentially adding sequences is continued until an intermediate tree is reconstructedcontaining all n taxa.
The puzzling step itself is repeated several thousand times with different input orders.
3.2.2.3. Consensus Step
From the collection of all intermediate trees a majority-rule consensus tree is constructed,
resolving all bipartitions which occur in more than 50% of the intermediate trees (formore details on consensus methods see 7.2.2.1 as well as McMorris and Neumann, 1983).
The percent occurrence of the bipartitions among the intermediate trees gives a support
value for the corresponding branch.
Finally, maximum likelihood branch lengths and the log-likelihood are computed forthe consensus tree following the approach of Felsenstein (1981) as described in 2.2.3.

4. Improvement in Complexity of the
Puzzling Step of QP
4.1. Introduction
Large growth rates of running times with increasing input size is a serious obstacle forthe application of complex methods to large datasets.
When running time gets large more efficient algorithms are necessary. The enormous improvement that can be achieved by efficient algorithms has been impressively
demonstrated in computer science textbooks, e.g., for sorting and searching algorithms(Cormen et al., 2001; Knuth, 1997; Ottmann and Widmayer, 2002).
The efficiency of computer algorithms can be measured by running benchmarks on acomputer to compare the performance of two algorithms. Benchmarks are a possible
and necessary way to judge the efficiency of, e.g., a parallelized algorithm compared tothe sequential version (cf. chapter 5). To measure the global efficiency of computer algorithms, benchmarks are by far too machine dependent. Hence, the general complexityof algorithms is scrutinized to reveal their efficiency. In case of algorithms with similaror equal complexity, benchmarks can help to find the preferable method.
In general two kinds of efficiency are investigated, the consumption of memory andruntime. Since one is interested in the consumption of resources for large input sizes,
complexity measures have been introduced to describe the growth of computational andmemory needs with respect to the input size.
This chapter is concerned with the complexity of the puzzling step introduced inchapter 3. For this purpose, the relevant complexity measures will be presented in
the next section. We then will examine the complexity of the original puzzling step.
Afterwards, three more efficient algorithms will be described, two novel ones and one
suggested by Vincent Berry (Strimmer, personal communication). Finally, the runtimeof all four algorithms will be compared and the difference will be discussed.

4.2. Complexity Measures
The complexity of an algorithm is usually described by asymptotic bounds using the O,Ω, and Θ-notation (Cormen et al., 2001, chap. 3). The three notations are defined asfollows.

25

26

4. Improvement in Complexity of the Puzzling Step of QP

a

c

b

c2 g(n)c g(n)
f(n)

f(n)

f(n)
c g(n)n0

n0

f (n) = O(g(n))

f (n) = Ω(g(n))

c1 g(n)n0

f (n) = Θ(g(n))

Figure 4.1.: The three complexity notations O, Ω, and Θ
The O-notation describes an asymptotic upper bound of an algorithm and is definedbyO(g(n)) = {f (n) : there exist c, n0 ≥ 0,such that 0 ≤ f (n) ≤ cg(n),for all n ≥ n0 }(4.1)

This equation means function f (n) belongs to the set O(g(n)), if there exists a positive
constant c, such that cg(n) is an upper bound for f (n) for all n ≥ n0 (Fig. 4.1a). f (n)might, for example, describe the memory consumption of an algorithm.
The Ω-notation defines an asymptotic lower bound.
Ω(g(n)) = {f (n) : there exist c, n0 ≥ 0,such that 0 ≤ cg(n) ≤ f (n),for all n ≥ n0 }(4.2)

If there is a positive constant c, such that cg(n) is a lower bound of f (n) for all n ≥ n0
(Fig. 4.1b), f (n) belongs to the set Ω(g(n)).
The Θ-notation is defined asΘ(g(n)) = {f (n) : there exist c1 , c2 , n0 ≥ 0,such that 0 ≤ c1 g(n) ≤ f (n) ≤ c2 g(n),for all n ≥ n0 }(4.3)

This means function f (n) belongs to the set Θ(g(n)), if there exist positive constants
c1 , c2 , n0 , such that f (n) can be ’sandwiched’ between c1 g(n) and c2 g(n) for all n ≥ n0
(Fig. 4.1c). Note, if f (n) belongs to Θ(g(n)), this implies that it also belongs to O(g(n))
as well as to Ω(g(n)).
O and Ω can be taken as the worst-case or best-case complexity, respectively.

4.3. Complexity of the Original Puzzling Step

27

Examples An algorithm that adds up all entries of an (n × n) matrix should be oforder Θ(n2 ), because every entry has to be touched exactly once.
2

Given a symmetric matrix, only n 2+n entries, that is, one triangular part, have to beaccessed. This, however, is still included in Θ(n2 ).
Searching whether a zero contained in such a matrix, in worst case, all entries haveto be touched once if there is no zero or one in the last entry. In the best case a zero isfound in the first entry. Hence, such an algorithm should be of order O(n2 ) and Ω(1).
Most often the upper bound, denoted by the O-notation, is used to characterize analgorithm. Important growth rates for algorithms with respect to input size n are
O(1)O(log n)O(n)O(n log n)O(n2 )O(n3 )O(nx )O(xn )

constant time (no growth dependent on input size n)logarithmic growthlinear growthn log n growthquadratic growth
cubic growthpolynomial growth
exponential growth.

One algorithm is considered more efficient than another if it has lower a growth rate.
For more details refer to computer science textbooks (e.g., Cormen et al., 2001; Ottmann
and Widmayer, 2002).

4.3. Complexity of the Original Puzzling Step
The original algorithm for the puzzling step (cf. 3.2.2) was implemented as described by
Strimmer and von Haeseler (1996). A pseudo-code listing of the algorithm (Listing 4.1)
will be used to examine its complexity. In the listing, the lines and loops important forthe measurement of the complexity are numbered (e.g., 4.1c).
Length of Penalty Paths and the Diameter of a Tree Central to the original algorithm is the incrementation of the edge penalties (step 4.1c) on the path between two
penalty neighbors a and b, the two leaves not being neighbors of the insertion taxon X
in a preferred quartet topology (cf. 3.2.2 and Fig. 3.5).
The complexity of this step equals the length of the path between a and b, because
direction vectors are used to look up the direction at any node to any other node in thetree. The entries of the vectors indicate directions to any internal or external node. The
directions to the internal nodes are necessary to be able updating the direction to newly
inserted edges easily. For the following it is assumed that tree Ti with i leaves has been
constructed so far. The maximal path length between penalty neighbors is bounded by
the diameter of the tree Ti , which is defined as the longest path between any two leavesof the tree.

28

4. Improvement in Complexity of the Puzzling Step of QP

Listing 4.1: Original Puzzling Step Algorithm: Set of taxa S = {s1 , s2 , . . . sn };Set of preferred quartet topologies for all quartets in S
Result : One binary tree Tn containing all n taxa in S
begin
Set the sequences in S into a random order s1 , s2 , . . . sn ;Set tree T4 one preferred topology of {s1 , s2 , s3 , s4 };Initialize path directions in T4 ;4.1afor i = 4 to (n − 1) do
X = si+1 ;
4.1bforeach quartet qX ∈ Q(X,Ti ) do
Let a and b be the penalty neighbors in the preferred topology of qX ;4.1c
Add a penalty of 1 to all edge penalties Φ(ε) for all edges ε on the penalty
path between a and b in Ti ;4.1d
Find the set of edges Emin (Ti ) with minimal penalty;
Insert X into one edge e ∈ Emin (Ti ) ;Call the new tree Ti+1 ;4.1eUpdate path directions in Ti+1 ;return Tn ;
end
Data

The longest possible diameter for a tree with i leaves exists in a linearized tree topology(cf. Fig. 2.4d). Its diameter is exactly i − 1, since all i − 3 internal edges are on the pathplus two external branches. Hence the upper bound for the diameter is O(i).
The shortest possible diameter of a binary tree exists in a completely balanced tree
(Fig. 2.4c). It is well-known (e.g., Ottmann and Widmayer, 2002) that the maximum
height hmax of a balanced binary rooted tree (as in Fig. 2.4a) with i leaves is hmax ≤dlog2 ie, where the height h is defined as the number of edges between a node and theroot. Hence, in this case the height h is bound by Θ(log i). Therefore, the diameter ofsuch rooted trees is twice its height, implying Θ(log i) for the diameter. The same holds
for the unrooted tree the diameter being one edge shorter because the central edge is
not divided by a root (Fig. 2.4c). Accordingly, the complexity of step 4.1c incrementing
the penalty path is bound by O(i).

Runtime complexity Step 4.1c is repeated for each of the 3i relevant quartets (loop
4.1b), leading to complexity 3i · O(i) ⊂ O(i4 ) for loop 4.1b. Finding the set of minimaledges Emin (Ti ) has complexity O(i), since all 2n − 3 edges have to be examined (step4.1d). After inserting X, the directions pointing to X and the new internal node have tobe added by inserting the direction to the insertion edge in constant time per edge (step
4.1e. The direction vector of one of the adjacent nodes is copied to the newly inserted

4.4. More Efficient Puzzling Step Algorithms

29

node (O(i)). The linear complexity of this update (step 4.1d) is superseded by the O(i4 )complexity of loop 4.1b. Since we add the n taxa sequentially (loop 4.1a), the overallruntime complexity of the original puzzling step algorithm is O(n5 ).

Memory complexity The memory needed for general purposes such as tree topologies
and leaf set permutations are of linear complexity. This remains the same in all following
algorithms and will, thus, not be considered any further. As mentioned above, direction
vectors are used at all internal nodes to efficiently determine paths through the tree.
The memory complexity of the direction vectors is Θ(n2 ).

4.4. More Efficient Puzzling Step Algorithms
A main improvement in complexity compared to the original algorithm can be achievedby avoiding the tree traversal for each quartet during the insertion step.
Here two new algorithms are presented which reduce the complexity of the puzzling
step. Both algorithms are based on a so-called penalty neighborhood matrix N = (Na,b )with a, b ∈ L(Ti ), the leafset
 of Ti . The entries of N count the number of times twoitaxa are neighbors in the 3 quartet tree topologies relevant to insert a next taxon X
into Ti
Na,b = |{ ab|cX | a, b, c ∈ L(Ti )}|(4.4)

The two algorithms differ in the way the penalty Φ(e) for an edge e is derived from thepenalty neighborhood matrix N .
Additionally, an algorithm suggested by Vincent Berry (Strimmer, personal communication)is presented that is based on the MRCA (most recent common ancestor) relationships of the penalty neighbors.

4.4.1. Split-Based Puzzling Step Algorithm
The following algorithm is sketched in Listing 4.2. To compute the penalty of an edge e
only the sum of penalty paths passing through e induced by any two penalty neighbors on
different sides of that edge have to be considered. The sum of all paths between any two
penalty neighbors are computed in the penalty neighborhood matrix N (Eq. 4.4). Each
incrementation has complexity O(1) which is performed for the 3i relevant quartets.
Accordingly, computing N (loop 4.2b) has complexity Θ(i3 ), since symmetrizing thematrix (loop 4.2c) needs Θ(i2 ) assignments, which does not increase the complexity. N
is symmetrized to avoid sorting every pair of leaves when accessing N .

30

4. Improvement in Complexity of the Puzzling Step of QP

Listing 4.2: Split-based Puzzling Step
: Set of taxa S = {s1 , s2 , . . . sn };Set of preferred quartet topologies for all quartets in S
Result : One binary tree Tn containing all n taxa in S
begin
Set the sequences in S into a random order s1 , s2 , . . . sn ;Set tree T4 one preferred topology of {s1 , s2 , s3 , s4 };Initialize split vectors in T4 ;4.2afor i = 4 to (n − 1) do
Clear (i × i) penalty neighborhood matrix N ;X = si+1 ;
4.2bforeach quartet qX ∈ Q(X,Ti ) do
Let a = sx and b = sy (x < y) be the penalty neighbors in the preferred
topology of qX ;Increment Na,b ;
Data

4.2c
4.2d
4.2e

foreach x < y ≤ i do Nsy ,sx = Nsx ,sy ;foreach edge e in Ti do
foreach leaf a ∈ e doforeach leaf b ∈ e do
Increment Φ(e) by Na,b ;4.2f

find the set of edges Emin (Ti ) with minimal penalty;
insert X into one edge e ∈ Emin (Ti );call the new tree Ti+1 ;
4.2h
update split vectors in Ti+1 ;return Tn ;
end
4.2g

The penalty Φ(e) of an edge e is collected (loops 4.2e and 4.2f) summing up the joint
penalties Na,b of each penalty path from a and b leading from any leaf a on the one sideof e to any leaf b on the other side of e,Φ(e) =XX

Na,b .

(4.5)

a∈e b∈e

The computation of one Φ(e) needs between i−1 for external edges and d(i/2)2 e additions, that means loops 4.2e and 4.2f together have complexity O(i2 ). The computationof the penalties for all 2i − 3 edges has overall complexity O(i3 ) (loop 4.2d).
Finding the set of minimal penalty edges Emin (Ti ) (step 4.2g) is again of order O(i),but can be collected while computing Φ(e).

4.4. More Efficient Puzzling Step Algorithms

31

To be able to look up the splits at each edges, we use two split arrays of joint size n foreach edge, which have to be updated once after an insertion (step 4.2h). After insertion
we have to add the newly inserted taxon to the existing arrays in constant time per edge
in a tree traversal, and all existing taxa have to be added to the cluster vectors at the
two newly created edges again in O(i).
Overall runtime complexity The highest complexity within the insertion loop 4.2a isO(i3 ) for the computation of all edge penalties as well as for the construction of thepenalty neighborhood matrix N , which supersede the other complexities. Since we add
all taxa one by one, the algorithm’s overall runtime complexity is O(n4 ).
Memory complexity In comparison to the original algorithm, the split-based algorithm
requires additional memory for the penalty neighborhood (n × n)-matrix N and cluster
vectors with a joint size of n entries at every edge. The data structures both havecomplexity Θ(n2 ). Hence, the memory complexity is equal for both algorithms. Note,the directional arrays of the original algorithm are not needed any more.

4.4.2. Recursive Puzzling Step Algorithm
As described above, the split-based algorithm reduces the complexity of the puzzling stepby one order of magnitude. Yet further improvement can be achieved. While examining
all 3i relevant quartets is necessary for the insertion, computing the edge penalties canbe solved more efficiently by using the following recursive approach (Listings 4.3 and
4.4).
The recursive algorithm again uses the penalty neighborhood matrix N which is constructed as before (Eq. 4.4). Additionally, the algorithm uses an external node penalty
vector E = (Ea ) containing the number of penalty paths starting at each leaf a. E iscalculated as follows.
Ea =

X

Na,b

a, b ∈ L(Ti )

(4.6)

b6=a

Starting from N and E, the edge penalties are computed via post-order traversal (calling
the recursive function, Listing 4.4, in step 4.3d) in the recursive puzzling step algorithm(Listing 4.3). The tree Ti is assumed to be rooted arbitrarily at the external branchleading to s1 , the first sequence in the random permutation. Note, that the placementof the root does not change the result. Yet the root is necessary to introduce a starting
point for the recursion.
Before analyzing the complexities of the algorithm sketched in Listings 4.3 and 4.4,
the recursive procedure will be described in detail using an example.

32

4. Improvement in Complexity of the Puzzling Step of QP

Table 4.1.: Quartets qH ∈ Q(H,T7 ) with L(T7 ) = {A, B, C, D, E, F, G} (cf. Fig. 4.1) andtheir supported topologies for the example.
{A, B, C, H}{A, B, D, H}{A, B, E, H}{A, B, F, H}{A, B, G, H}{A, C, D, H}{A, C, E, H}{A, C, F, H}{A, C, G, H}

AC|BH ∗
AB|DH
AB|EH
AB|F H
AB|GH
AC|DH
AC|EH
AC|F H
AC|GH{A, D, E, H}{A, D, F, H}{A, D, G, H}{A, E, F, H}{A, E, G, H}{A, F, G, H}{B, C, D, H}{B, C, E, H}{B, C, F, H}

rootedge ee(1,A)
e(1,B)e(root,C)e(3,D)e(3,E)
e(4,F )e(4,G)e(2,1)e(root,2)e(5,3)e(5,4)e(2,5)
root

13
2
12
17

5

13

H
11

1

C

1

3

4

12

13

10

11

6

5

A

B

D

E

F

G

DE|AH
DF |AH ∗AD|GH
AE|F HAE|GH
F G|AH
BC|DH
BC|EH
BC|F H{B, C, G, H}{B, D, E, H}{B, D, F, H}{B, D, G, H}{B, E, F, H}{B, E, G, H}{B, F, G, H}{C, D, E, H}{C, D, F, H}

BC|GH
DE|BH
BD|F H
BD|GH
BE|F H
BE|GH
F G|BH
DE|CH
CD|F H{C, D, G, H}{C, E, F, H}{C, E, G, H}{C, F, G, H}{D, E, F, H}{D, E, G, H}{D, F, G, H}{E, F, G, H}

CD|GH
CE|F H
CE|GH
F G|CH
DE|F H
DE|GH
F G|DHF G|EH

CLP (C e (e)) SMS (C e (e)) Φ(e) cluster C e (e)
12
0
12
{A}
13
0
12
{B}
13
0
13
{C}
10
0
10
{D}
11
0
11
{E}
6
0
6
{F }
5
0
5
{G}
25
8
17
{A, B}
57
44
12{A, B, D, E, F, G}
21
10
11{D, E}
11
10
1{F, G}
32
20
12{D, E, F, G}
7070
0{A, B, C, D, E, F, G}

Figure 4.2.: Edge penalties and penalty sums SMS and CLP computed recursively on
the tree T7 with N and E (Eq. 4.7) derived from the quartet topologies inTab. 4.1.
4.4.2.1. Example
One insertion step of the recursive Puzzling Step algorithm will be exemplarily performedon the tree T7 in Fig. 4.2 using the quartet topologies from Tab. 4.1.
According to Eqs. 4.4 and 4.6 the penalty neighborhood matrix N and its row sums,i.e., the external node penalty vector E are computed.






N =




A
4
5
1
2
0
0

B
4
4
2
2
1
0

C
5
4
2
2
0
0

D
1
2
2
5
0
0

E
2
2
2
5
0
0

F
0
1
0
0
0
5

G
0
0
0
0
0
5
-











A
B
C
D
E
F
G




12
 13 

 13 
E = 10  11 
 6 
5

A
B
C
D
E
F
G

(4.7)

4.4. More Efficient Puzzling Step Algorithms

33

Using Ea , which contains all penalties induced by penalty paths of leaf a, the penaltyof the external edge e(•,a) of leaf a can be setΦ(e(•,a) ) = Ea .

(4.8)

Climbing to an inner edge e(•,v) that descends to an inner node v, the edge penalty
Φ(e(•,v) ) can be computed as follows. All penalty paths originating at any leaf of theedge’s cluster C e (e(•,v) ) are added. The sum is called cluster leaf penalty CLP .
X

CLP (C e (e(•,v) )) =

Ea

(4.9)

a∈C e (e(•,v) )Note, that all penalty paths that start and end within the cluster are counted twice. To
obtain the true edge penalty we have to subtract all penalty paths remaining inside thecluster twice. This number is given by the submatrix sum of N
XX

k1 ∈C e (e(•,v) )k2 ∈C e (e(•,v) )SMS (C e (e(•,v) )) =Nk1 ,k2 .

(4.10)

Consider edge e(2,5) from Fig. 4.2. The edge’s cluster consists of the union of thedescending edges’ clusters, i.e., C e (e(2,5) ) = C e (e(5,3) ) ∪ C e (e(5,4) ) (cf. table in Fig. 4.2).
Hence, CLP (C e (e(2,5) )) is efficiently computed by adding the CLP s of the descending
branches, avoiding double calculations. Also SMS (C e (e(2,5) )) is computed more efficiently than summing up the whole submatrix of N colored grey in Eq. 4.7. The lightgrey parts have already been computed for the descendant edges. Thus, only one of thedark grey submatrices has to be summed up, because of N being symmetric:
SMS (C e (e(2,5) )) = SMS (C e (e(5,3) )) + SMS (C e (e(5,4) ))
X
X
+ 2·
Na,b
a∈C e (e(5,3) )(4.11)b∈C e (e(5,4) )Finally, by subtracting SMS from CLP the edge penalty
Φ(e(2,5) ) = CLP (C e (e(2,5) )) − SMS (C e (e(2,5) )).
of edge e(2,5) is computed.

(4.12)

34

4. Improvement in Complexity of the Puzzling Step of QP

Listing 4.3: Recursive Puzzling Step
: Set of taxa S = {s1 , s2 , . . . sn };Set of preferred quartet topologies for all quartets in S
Result : One binary tree Tn containing all n taxa in S
begin
Set the sequences in S into a random order s1 , s2 , . . . sn ;Set tree T4 one preferred topology of {s1 , s2 , s3 , s4 };Root Ti−1 arbitrarily at the external edge of s1 , then called eroot ;Initialize cluster vectors in T4 ;4.3afor i = 4 to (n − 1) do
Clear (i × i) penalty neighborhood matrix N ;Clear external node penalty i-vector E;
X = si+1 ;
4.3bforeach quartet qX ∈ Q(X,Ti ) do
Let a = sx and b = sy (x < y) be the penalty neighbors in the preferred
topology of qX ;Increment Na,b ;
Data

4.3c

foreach x < y ≤ i do
Nsx ,sy = Nsy ,sx ;Increment Esy by Nsy ,sx ;Increment Esx by Nsy ,sx ;Φ(eroot ) = SubtreePenalty(eroot , T , N , E);
if Φ(eroot ) 6= Es1 then error
4.3eFind the set of edges Emin (Ti ) with minimal penalty;
Insert X into one edge e ∈ Emin (Ti ) with minimal penalty;
Call the new tree Ti+1 ;4.3f
Update cluster vectors in Ti+1 ;return Tn ;
end

4.3d

4.4.2.2. Complexity
Runtime Complexity The recursive algorithm again uses a neighborhood matrix N
which is built and symmetrized the same way as before in Θ(i3 ) (loop 4.3b). Symmetrizing N and computing the row sums E (loop 4.3c) has complexity Θ(i2 ). External edge
penalties Φ(e(•,a) ) = Ea are assigned in constant time for each of the i external edges(condition 4.4a).
Ascending to an internal edge ε the penalty Φ(ε) is computed by adding the sums of
the penalty paths CLP of the clusters C e (e(•,l) ) and C e (e(•,r) ) of the descending edges e(•,l)
and e(•,r) (step 4.4e) in O(1). From this sum the penalties (steps 4.4b to 4.4d) inducedby paths among the members of the cluster C e (ε) of edge ε has to be subtracted (step4.4. More Efficient Puzzling Step Algorithms

35

Listing 4.4: SubtreePenalty(e(v,w) , T , N , E)
Data

: Current edge e(v,w) leading away from root, ending at node w;
Binary tree topology T ;Penalty Neighborhood Matrix N ;External Node Penalty Vector E
Result : Edge penalty Φ(e(v,w) ) of edge e(v,w) ;Cluster C n (w);
Sum of paths CLP (C n (w)) starting at leaves in cluster C n (w);
Sum of paths SMS (C n (w)) within cluster C n (w)
begin

4.4a

4.4b
4.4c

if w = external node thenΦ(e(v,w) ) = Ew ;CLP ({w}) = Ew ;SMS ({w}) = 0;
return Φ(e(v,w) ), C n (w), CLP (C n (w)), SMS (C n (w));
else
l = left descendant of w;(Φ(e(w,l) ), CLP (C n (l)), SMS (C n (l))) =SubtreePenalty(e(w,l) , T , N , E);r = right descendant of w;(Φ(e(w,r) ), CLP (C n (r)), SMS (C n (r))) = SubtreePenalty(e(w,r) , T , N ,E);
SMS (C n (w)) = 0;
foreach a ∈ C n (r) do
foreach b ∈ C n (l) do
Increase SMS (C n (w)) by Na,b ;SMS (C n (w)) = 2 · SMS (C n (w)) + SMS (C n (l)) + SMS (C n (r));
CLP (C n (w)) = CLP (C n (l)) + CLP (C n (r));Φ(e(v,w) ) = CLP (C n (w)) − SMS (C n (w));
return Φ(e(v,w) ), CLP (C n (w)), SMS (C n (w));

4.4d
4.4e
4.4f

end

4.4f), an operation with complexity O(i2 ) executed once per edge. Although this step
still has complexity O(i2 ) it reduces the runtime as described above. Finding the setof minimal penalty edges Emin (Ti ) (step 4.2g) costs again O(i), but can be collected
immediately while computing Φ(e).
After inserting the next leaf, the newly inserted inner edge has to get an initial clusterby copying it from the insertion edge in O(i) (step 4.3f). In addition, the inserted leaf
has to be added to the clusters of the edges on the path from the insertion edge to theroot in O(i) (step 4.3f).

36

4. Improvement in Complexity of the Puzzling Step of QP

The most expensive parts in the insertion loop 4.3a have complexity O(i3 ) per insertion, namely, the construction of N as well as computing the submatrix sums SMS
recursively for all edges (SubtreePenalty, Listing 4.4, called in line 4.3d). Since we
add all taxa one after the other, the overall complexity of the recursive algorithm isO(n4 ).
Memory complexity In comparison to the original algorithm, the recursive algorithm
requires additional memory for the penalty neighborhood (n × n)-matrix N , the external
node penalty n-vector E, and a vector of n entries per edge to store the cluster of theedge. Hence, the recursive algorithm has memory complexity of Θ(n2 ) as well.

4.4.3. Berry’s MRCA-Based Puzzling Step Algorithm
Another O(n4 ) algorithm has been proposed by Vincent Berry (Strimmer, personal communication). Berry’s algorithm calculates the edge penalties using the MRCA relationsof the penalty neighbors and, hence, is referred to as MRCA-based puzzling step algorithm. In this algorithm, all penalty paths between two leaves a and b are collected in
an n-vector at their MRCA. The entries of the vectors are then propagated down to therespective leaves in a pre-order traversal of the rooted tree (calling the recursive function
Listing 4.6 in step 4.5f).
The Algorithm As in the recursive algorithm the tree Ti is rooted arbitrarily at theexternal branch leading to s1 . Again, the placement of the root does not change theresult, but it is needed to introduce directions to define the MRCAs. In the beginning,
the first tree and the MRCA matrix
is set up in constant time. To add a sequence X toia previously built tree Ti , all 3 relevant quartets containing X and three leaves alreadyin Ti are examined (loop 4.5c). Let a and b be the penalty neighbors of such a quartet,then the entries of the penalty propagation vector V at the inner node MRCA(a, b) for
the nodes a and b are incremented (loop 4.5c). Afterwards, the entries of the penalty
propagation vectors an inner node are propagated down to the corresponding leaves ina pre-order traversal (called in step 4.5f). All entries of the propagation vector V areadded to the entries at the descendant node’s vector V , if the corresponding leaf belongs
to the cluster of this descendant node (steps 4.5d, 4.6a, and 4.6c). The penalty Φ(e) ofan edge e is computed by summing up all entries of V attached to the parental node ofe which correspond to a leaf in the cluster C e (e) of e (steps 4.5e, 4.6b, and 4.6d). The
taxon X is added at an edge with minimum penalty. (step 4.5g).
Runtime Complexity To initially fill the penalty propagation vectors requires evaluating all relevant quartets and, thus, has complexity O(i3 ) (loop 4.5c). Propagating
penalties from one node to the other (steps 4.5d, 4.6a, and 4.6c) has complexity O(i) as
has the computation of the penalty of an edge (steps 4.5e, 4.6b, and 4.6d), due to thenumber of possible leaves in a cluster. Since the propagation is performed recursively on4.4. More Efficient Puzzling Step Algorithms

37

Listing 4.5: MRCA-based Puzzling Step
: Set of taxa S = {s1 , s2 , . . . sn };Set of preferred quartet topologies for all quartets in S
Result : One binary tree Tn containing all n taxa in S
begin
Set the sequences in S into a random order s1 , s2 , . . . sn ;Set tree T4 one preferred topology of {s1 , s2 , s3 , s4 };Root T4 arbitrarily in the external edge to s1 ;Initialize MRCA (n × n)-matrix M for the rooted T4 ;4.5afor i = 4 to (n − 1) do
Data

4.5b

4.5c

4.5d

4.5e

Clear MRCA penalty (i−3 × i−3)-matrix V ;X = si+1 ;foreach quartet qX ∈ Q(X,Ti ) do
Let sx and sy (x < y) be the penalty neighbors in the preferred topologyof qX ;
Increment Vsx ,MRCA(sx ,sy ) ;Increment Vsy ,MRCA(sx ,sy ) ;l = succ(root)\s1 ;foreach d ∈ C n (l) do
Increase Vd,l by Vd,root ;Increase Φ(e(•,l) ) by Vd,l ;Φ = PropagatePenalties(e(•,l) , V , Φ);
4.5g
Find the set of edges Emin (Ti ) with minimal penalty;
Add X to one edge e ∈ Emin (Ti ) with minimal penalty;
Call the new tree Ti+1 ;4.5h
Update MRCA matrix for X with all leaves in L(Ti );4.5i
Update cluster vectors in Ti+1 ;return Tn ;
end
4.5f

all edges, the whole calculation of the 2i − 3 edge penalties has complexity O(i2 ). The
set of minimal edges (loop 4.5g) is collected during the traversal in O(1) at each edge.
After inserting X into Ti , X has to be added to the clusters of the nodes (step 4.5i)
between the newly inserted inner node and the root (complexity O(i)). Also, X has
to be added to the MRCA matrix (step 4.5h). The update follows the nodes between
the newly inserted taxon and the root, where the current node is the MRCA of X andany leaf in the cluster of the current edge which was not in the cluster of the previous
inner node. Since the clusters are known for the right and left descendant of any innernode and only i entries have to be added, the update has complexity O(i). The largest
complexity during one insertion step is O(i3 ) examining all relevant quartets (loop 4.5c),

38

4. Improvement in Complexity of the Puzzling Step of QP

Listing 4.6: PropagatePenalties(e(•,w) , V , Φ)
Data

: Current edge e(•,w) leading away from root, ending at node w;
Tree topology T ;MRCA penalty matrix V
Result : Edge penalty Φ(e(•,w) ) of edge e(•,w)
begin
l = left descendant of w;
4.6aforeach a ∈ C n (l) do
Increase Va,l by Va,w ;4.6b
Increase Φ(e(w,l) ) by Va,l ;
4.6c

4.6d

r = right descendant of w;
foreach b ∈ C n (r) do
Increase Vb,r by Vb,w ;Increase Φ(e(w,r) ) by Vb,r ;if r internal node then PropagatePenalties(e(w,l) , V , Φ);if l internal node then PropagatePenalties(e(w,r) , V , Φ);
endwhich makes the overall runtime (for loop 4.5a) complexity O(n4 ).
Memory complexity The data structures maintained by the algorithm are the MRCA
matrix with Θ(n2 ) and the propagation vectors V at all inner nodes with Θ(n2 ). Hence,
the memory complexity of this algorithm is again Θ(n2 ) as in all other puzzling stepalgorithms proposed in this chapter.

4.5. Runtime Analysis
All improved algorithms presented here have the same reduced runtime complexity ofO(n4 ) while keeping the memory complexity at Θ(n2 ). To measure the performance
improvement of the different algorithms, benchmark tests have been performed for theoriginal as well as the improved algorithms.

4.5.1. Datasets
To test the runtime performance of the implementation of the four different puzzling
step algorithms, we used a variety of datasets with different numbers of taxa.
Datasets of different sizes were used: two biological amino acid datasets of elongation
factors, the first consisting of size 24 (EF-24 Schmidt and von Haeseler, 2003) with 915

4.5. Runtime Analysis

39

Table 4.2.: Results of the runtime analysis of the four different puzzling step algorithms:
the original algorithm, the MRCA-based, split-based, and recursive modifications. The best runtime is marked bold-face.
datasetseqs. length original MRCA-based split-based
recursive
tsplittMRCA
torig tMRCA
tsplit torig
trecur ttrecur
torig
orig
EF-24
24 915aa
32s
20s
0.63
20s 0.63
19s 0.59
EF-50
50 1117aa
1159s
610s
0.53 589s 0.51557s 0.48
GST-5050 1065aa
1191s
582s
0.49 567s 0.48537s 0.45
Sim-32-500
32
500nt
154s
91s
0.59
89s 0.58
84s 0.55
Sim-32-1000
32 1000nt
154s
91s
0.59
89s 0.58
83s 0.54
Sim-48-500
48
500nt
894s
477s
0.53 460s 0.51438s 0.49
48 1000nt
887s
472s
0.53 455s 0.51432s 0.49
Sim-48-1000
64
500nt
3146s 1614s
0.51 1541s 0.49 1475s 0.47
Sim-64-500
Sim-64-1000
64 1000nt
3100s 1548s
0.50 1480s 0.48 1415s 0.46

aligned columns and the other comprising 1117 aligned sites from 50 sequences (EF50). As a third biological dataset cluster 63525 from the SYSTERS database (Version
3; Krause et al., 2002, http://systers.molgen.mpg.de/) was used, consisting of 50
glutathione S-transferase proteins with 1065 aligned columns (GST-50).
Additionally, simulated datasets with different numbers of taxa (32, 48, 64) as well as
different sequence lengths (500nt, 1000nt) were used (for more details on the simulated
datasets refer to 5.3).

4.5.2. Benchmark Setup
To get even conditions, the initial random number seeds in any run were fixed, suchthat the runs only differed in the construction of the intermediate trees by the algorithms, but not in the choices made concerning quartet topologies or insertion edges. In
each benchmark run 10,000 intermediate trees were computed. The benchmarks were
performed several times and the runtimes were averaged and rounded to the next fullsecond.
The benchmarks were computed on a LINUX workstation with an AMD Athlon XP
2100+ CPU and 512MB memory.
Since storing the splits is done at the end of each single puzzling step, a portion of theruntime is unaffected by the puzzling step algorithms. It was not possible to determine
the exact split maintenance times, since – especially for small tree sizes – these are nearor even below the resolution of the hardware clocks used for the measurements.

40

4. Improvement in Complexity of the Puzzling Step of QP

4.5.3. Results
The results of the benchmark tests (Tab. 4.2) show a clear reduction of the runtime by
the three O(n4 ) algorithms of about 50% compared to the runtime of the original O(n5 )algorithm. The best reduction was achieved by the recursive algorithm with 59-45% of
the original algorithm’s runtime. It is followed by the split-based algorithm with 6348% runtime and then by Berry’s MRCA-based algorithm with 63-49%. Although, thepercentages are sometimes close to each other, the algorithms end up clearly ordered bytheir runtimes, namely the MRCA-based algorithm, the split-based algorithm, and therecursive algorithm, ordered by decreasing runtime for all benchmark datasets.

4.6. Discussion
The results in Tab. 4.2 convincingly demonstrate the positive effect of using more efficient algorithms. As described in 4.5.3 all O(n4 )-algorithms show reduced runtimesof about 50% compared to the original algorithm for all benchmark datasets. Furthermore, the three improved algorithms show a clear order of their runtimes. Although
the implementation of the MRCA-based algorithm substantially reduced the runtime it
was outperformed by the split-based as well as the recursive puzzling step algorithm. As
expected, the split-based algorithm show longer runtime compared to the recursive implementation, which is due to the computational overhead calculating the edge penalties
from the penalty neighborhood matrix N (cf. 4.4.2).
In particular, if we consider the operations add, assign, multiply, and subtract, computing the calculation of SMS (C e (e(2,5) )) (Eq. 4.11) in the example on page 32 needs
four fix operations and four variable operations within the sums (variable means thatthe number depends on the cluster sizes). The calculation of CLP (C e (e(2,5) )) needs two
fix operations. Two more operations are needed to compute Φ(e(2,5) ) (Eq. 4.12). In theend eight fix operations and four variable ones are needed, compared to twelve variable
operations in the split-based algorithm. Note, that variable operations take longer time
due to the counters and comparisons within the addition loops. Hence, although thenumbers of operations counted are equal in this case a speedup can be gained. Additionally, a lot of dereferencing is avoided. The difference increases dramatically in larger
trees. For instance, the penalty of an edge joining two clusters of size 5 and 6 in a treewith 30 leaves, needs 30 variable and again 8 fix operations, compared to 19 · 11 = 209
variable operations in the split-based algorithm.
Furthermore, the results from the simulations also show that the sequence lengths have
only marginal effect on the runtimes measured. This is obvious, since the puzzling stepis based on the ML quartet topologies. The fluctuation between the 500nt and 1000nt
datasets are mainly to be attributed to a lower degree of information in the shorter
sequences. Subsequently, more different intermediate trees are reconstructed from theless resolved quartets. This increases the number of occurring splits which leads to aslightly higher overhead maintaining the increased number of splits.

4.6. Discussion

41

The results also show that the percentage of the remaining runtime decreases for
growing numbers of taxa, which is due to the difference of one order of magnitude
between the still polynomial complexities. The improved runtime complexity decreases
the runtime of all datasets analyzed, and thus, enables the analysis of large datasets inless time.

42

4. Improvement in Complexity of the Puzzling Step of QP

5. Parallelized Quartet Puzzling
5.1. Introduction
The often exponential complexities of the applications lead to the use of heuristics.
However, even heuristics sometimes do not succeed in reducing the runtime sufficiently
to enable the analysis of large datasets within reasonable time. Therefore, parallel
computing techniques emerged as valuable tools to reduce the runtime of applications.
Parallel computing has first been introduced to pairwise alignment algorithms forsequence comparison (Edmiston and Wagner, 1987; Lander et al., 1988; Huang, 1989).
Since the early 1990s, it has been applied to phylogenetic analysis (Hagstrom et al., 1992;
Olsen et al., 1994, fastDNAml) by parallelizing the DNAML program (Felsenstein, 1981).
Until today parallel processing has entered most areas of modern sequences analysis
like database searching (Jülich, 1995; Pedretti et al., 1999), RNA structure prediction(Nakaya et al., 1995; Shapiro et al., 2001), multiple sequence alignment (Kleinjung et al.,2002), and recently also contig assembly (Ness et al., 2002), and gene selection for tissueclassification (Liu et al., 2001).
Quite a variety of parallel platforms and techniques exist. Among the first parallel
computers were vector computers which contain an array of processing elements (PEs)that are able to simultaneously execute the same instruction (SIMD computers – single
instruction multiple data). More abundant nowadays are MIMD computers (multiple
instructions multiple data) that are able to compute different tasks simultaneously on
the different PEs/CPUs. The main types of MIMD platforms are multi-processor sharedmemory computers, also known as SMP (scalable multi-processor), distributed memorycomputers, and clusters of workstations (COWs). While all CPUs of a shared-memorycomputer have access to a global (shared) memory, CPUs in distributed memory computers all have their own ’private’ memory and communicate via a fast interconnection
device. A COW comprises a number of workstations via a local network (e.g., Ethernet). Employing workstations with multiple processors, characteristics of both SMP anddistributed memory platforms are combined.
Two main parallel programming paradigms exist, namely, shared memory and message passing. Shared-memory programming makes use of the shared memory either byaccessing the same copy of data and variable with all processes or by running several
subprocesses (threads) on the same program code in memory or both. In message passing, the processes communicate via explicit messages over the interconnection facilities.
While the former is clearly restricted to shared memory architectures, the latter can be

43

44

5. Parallelized Quartet Puzzling

implemented on all platforms mentioned. Thus, message passing gains high flexibility
and portability. For more details see the review by Trelles (2001) on the use of parallel
computing in bioinformatics and parallel platforms.
Although a number of parallelized phylogenetic programs have been published to date,mainly the DNAML program (Felsenstein, 1981) for the ML based analysis of DNA data
was subject to parallelization attempts (Olsen et al., 1994; Schmidt, 1996; Trelles et al.,
1998; Stewart et al., 2001; Stamatakis et al., 2002).
This chapter presents the efficient parallelization of the Quartet Puzzling algorithm
implemented in TREE-PUZZLE (cf. chapter 3) employing Message Passing Interface(MPI, Snir et al., 1998; Gropp et al., 1998), a highly portable de facto standard inparallel programming. In the next section, the default parameter estimation of TREEPUZZLE will be explained. The runtime profile of the sequential program will be analyzed, and the resulting parallelization scheme will be presented. The scalability of theparallel program will be tested, and the results will be discussed with respect to otherparallelized ML-based tree reconstruction methods. Finally, a short outlook on further
enhancements will be given.

5.2. Parameter Estimation for Evolutionary Models
As explained in 2.2.3, ML tree reconstruction methods require a model of the evolutionary process. The parameters that specify these substitution models need to be estimatedfrom the data sets.
To account for the fact that positions in an alignment may evolve with different rates,
each site may also get its own rate specific factor (cf. 2.1.4). Typically one assumes that
rates are distributed according to a Γ-distribution (Gu et al., 1995; Yang and Kumar,1994), where one parameter α describes the amount of rate heterogeneity. Also, theshape parameter α must be estimated from the data.
The collection of parameters can either be specified by the user or can be estimateddirectly from a multiple sequence alignment. The parameter estimation is one of thechallenging problems in maximum likelihood tree reconstruction. As for phylogenetic
analysis in general, it is assumed that incorporating more data increases the accuracyof the estimates. In TREE-PUZZLE the following heuristic estimation algorithm isimplemented (Strimmer, 1997).

5.2.1. Algorithm: Estimating Model Parameters
A single step of the heuristic estimation algorithm is described in Listing 5.1. Withan initial or already optimized parameter set M a distance matrix D and from D aneighbor-joining tree Tnj (Saitou and Nei, 1987) with branch lengths is inferred. WithTnj a new estimate of M is computed. The ML estimation is performed as explained in5.3. Datasets

45

Listing 5.1: Estimating Model Parameters in TREE-PUZZLE
Data: S = a set of n aligned sequences
Result : Model parameters M
begin
Assign biologically reasonable initial values to M ;repeatEstimate the ML distance matrix D for all n(n−1)
sequence pairs as described2by Adachi and Hasegawa (1996b);5.1bConstruct a neighbor-joining tree Tnj from D (Saitou and Nei, 1987);
5.1c
Fit the distances in D to Tnj by optimizing its branch lengths with the leastsquares method;
5.1d
On the tree with fixed branch lengths compute the maximum likelihood
estimate of parameters M ;until estimate of M does not change any more;
return M ;
end
5.1a

2.2.3, but assuming the tree Tnj fixed and optimizing the parameters of M to maximize
the likelihood. This procedure is repeated until the estimate of M does not change any
more.
The single estimation step as described in Listing 5.1 is repeated alternately for theparameters of the evolutionary model and rate heterogeneity parameters until the parameters set M converges or a maximal number of repetitions has been performed.

5.3. Datasets
For the analytic tree reconstructions as well as the benchmark tests we use simulated
datasets as well as a biological one.
We generated balanced phylogenetic trees with 32, 48, and 64 taxa. For these treesalignments of length 500 and 1000 nucleotides (nt) where generated using Seq-Gen (version 1.2.6 Rambaut and Grassly, 1997) and assuming the HKY model (Hasegawa et al.,1985) and Γ-distributed rates with α = 1.0 introducing intermediate heterogeneity
among sites (Gu et al., 1995). To make the results of the datasets comparable, thetree lengths from root to the leaves were set to be 0.2 substitutions per site.
One might argue that using balanced trees influences the runtime behavior of thealgorithm. Certainly, the runtimes of the parts differ relative to each other, but using
other tree topologies shows similar behavior. Hence, balanced trees are chosen becauseof their simple scalability.

46

5. Parallelized Quartet Puzzling

Table 5.1.: Runtime profile in percent of total runtime for ML parameter estimation.
Taxa
n
32
32
48
48
64
64
50

length500 nt1000 nt
500 nt1000 nt
500 nt1000 nt
289 aadistance matrix
mldistance
90.4 %
90.7 %
88.7 %
90.9 %
85.1 %
88.5 %
85.1 %

least square fit
lslength
2.6 %
1.4 %
6.6 %
3.5 %
11.3 %
7.4 %> 0.1 %In addition to the simulated data, 50 glutathione transferase sequences of the cluster
63,525 from the SYSTERS database (Version 3; Krause et al., 2002, http://systers.
molgen.mpg.de/) were used as a biological dataset. The alignment of the glutathionetransferase proteins has a length of 289 amino acids (aa).

5.4. Runtime Analysis of the Sequential TREE-PUZZLE
Runtime profiling (Tabs. 5.1 and 5.2) has been performed on a variety of datasets,
simulated as well as biological data to measure the time requirement of the different
parts in the tree reconstruction. This is done to locate the parts which demand large
portions of the overall runtime and, thus, are a prominent target for parallelization. Herethe results are presented for those datasets which are used later for the benchmarks (fordetails, cf. section 5.7).

5.4.1. Runtime of the Parameter Estimation
Runtime profiling (Tab. 5.1) shows that the repeated computation of the n(n−1)
pair2wise distances in step (5.1a) accumulates about 85-91% of the total computation time
during parameter estimation. The second largest part with up to 11.3% calculates theleast square branch lengths (step 5.1c) which performs (2n − 3)(n − 1) independent
computations per call.

5.4.2. Runtime of the Quartet Puzzling Algorithm
In the QP algorithm, the ML step and the puzzling step are by far most time-consuming(Tab.
 5.2). Both steps consist of many independent tasks, the computation of trees fornquartets and up to 50,000 computations of the intermediate trees. This makes both
4
steps a prominent target for parallelization.

5.5. Parallelizing TREE-PUZZLE

47

Table 5.2.: Runtime measurements for the different components of TREE-PUZZLE.
Taxa
length parameter number of ML stepnestimation
quartets
32
500 nt6.87 %
35,960 51.40 %321000 nt
8.45 %
35,960 65.14 %48
500 nt3.05 %
191,878 52.29 %48
1000 nt
3.61 %
191,878 69.44 %64
500 nt0.96 %
628,649 29.83 %641000 nt
1.27 %
628,649 43.15 %50
289 aa1.32 %
230,300 84.55 %# puzzling puzzling
steps
step
10,000 41.05 %10,000 25.84 %
10,000 44.43 %10,000 26.72 %25,000 68.99 %25,000 55.53 %10,000 13.51 %consensus
step0.68 %
0.56 %
0.21 %
0.23 %
0.21 %
0.05 %
0.62 %

5.5. Parallelizing TREE-PUZZLEFor portability of the code as well as applicability to different parallel platforms including
clusters of workstations (COWs), which are most common in molecular biology (Trelles,
2001), we decided to apply the message passing paradigm using the Message Passing
Interface standard (MPI, Snir et al., 1998; Gropp et al., 1998). Master/worker setups
in combination with dynamic scheduling algorithms are applied to gain flexibility andindependence from the underlying parallel platform. The choice of the applied schedulingalgorithms is discussed in section 5.6.
Due to serious differences in granularity (= the size of the independently computable
units in the sequential implementation) between the parameter estimation and the Quartet Puzzling algorithm both parts are optimized separately.

5.5.1. Parallelizing the Parameter Estimation
To reduce computation time, we applied a master/worker setup for both steps (Fig. 5.1).
The master process constructs the neighbor-joining trees (5.1b) and re-estimates theparameters (5.1d), while the computations of steps (5.1a) and (5.1c) are distributed
among the worker processes.

5.5.2. Parallelizing Quartet Puzzling
Again a master/worker scheme is used. The master process keeps track of the quartets tobe evaluated during the ML step and dispatches batches of tasks to the worker processes(Figure 5.1). To allow efficient transfer of the reconstructed quartet tree, the atomic
batch size has to be a multiple of two since the preferred topologies of one quartet isstored in 4 bit.
In the puzzling step the master process tells the workers for how many permutationsof taxa an intermediate tree needs to be reconstructed (Figure 5.1). The random permutations are based on SPRNG (http://sprng.cs.fsu.edu/), a parallel random number

48

5. Parallelized Quartet Puzzling
Multiple Sequence Alignment

ML DistLS LenScheduler

010

Evaluation of
all quartet trees

Puzzling Step:
Generating many
puzzling trees
from different
input orders

Scheduler

Scheduler

ML DistLS LenConsensus Step:
Strict majority rule consensus
from all puzzling trees
Evaluating branch lengths
and ML value of the new tree

010

Worker Process p

Master Process

ML Step:

Worker Process 1

Estimating
and sending
parameters

Maximum Likelihood Tree

Figure 5.1.: Parallelized workflow of the TREE-PUZZLE program. The upper part ofthe master process controls the initialization and the parameter estimation
where the least square branch length fitting (LS len) and ML distance estimation (ML dist) is performed by the worker processes. The lower partpresents the parallel Quartet Puzzling algorithm.
generator, to avoid dependencies of the workers’ generators. A list of all splits of theintermediate trees together with their absolute frequencies is returned to the master.
The number of splits occurring during the puzzling steps depends on the quality of thedata. Datasets with little phylogenetic information produce a higher number of different
intermediate trees and, hence, a larger number of splits.

5.6. Scheduling Algorithms
To ensure an even work load of the worker processes, we applied dynamic schedulingalgorithms (El-Rewini et al., 1994). Dynamic scheduling algorithms are important on

5.6. Scheduling Algorithms

49

heterogeneous COWs to produce a dynamically balanced work load. The effect of different scheduling algorithms on communication overhead and distribution of work loadon the p worker processors is investigated on a heterogeneous COW comprising SUN
workstations ranging from SPARC 20 to ULTRA 5.
For the parallelization scheme displayed in Figure 5.1 we tested algorithms which neither require knowledge about the network or allocation delays nor estimates of expected
runtime. The following algorithms were evaluated. The notation follows Hagerup (1997)
where the current status of the scheduler is denoted by Ξ and the return value is thesize of the next batch to schedule:• Self-Scheduling: SS(Ξ) = 1• Static Chunking: SC(Ξ) = NTNP
NR
• Guided Self-Scheduling: GSS(Ξ) = NP
NT
• Trapezoid Self-Scheduling: T SS(Ξ), T SS starts with batch size 2NP
decreasing
linearly to a final value of 1,

where NT is the total number of tasks, NR is the number of tasks not yet scheduled,and p denotes the number of worker processes.
To test the performance of the scheduling algorithms, we ran the TREE-PUZZLE
program several times with different scheduling algorithms on the heterogeneous SUN
COW and analyzed the according GANTT-chart diagrams produced by the program.
Since the presentation of the GANTT-charts is not possible we summarize the results.
Although the ’naive’ approach, self scheduling (SS), which sends all tasks one-by-one,
is able to produce a homogeneous work load, the communication and the subsequent
overhead are very high which leads to delays when sending task batches and results.
Static chunking (SC) performs nicely if all processors work at the same speed and if
the runtime of the tasks does not vary. In practice this is unrealistic for most parts of
TREE-PUZZLE. SC has the advantage of minimizing communication, the batch of aworker can even easily be pre-computed. On heterogeneous COWs, however, the slowest
processor becomes a bottleneck for the runtime of the whole set of tasks. Due to thevery small computation times per call in the least square branch length inference ofstep (5.1c), the communication overhead has to be kept as small as possible. Hence,, i.e. static chunking, are used. Only the initial dataprecomputed batch sizes of NT
NP
values are sent to the workers, and the results are then sent to the master process.
In preliminary tests, best performing were such scheduling algorithms which start witha large batch size, decreasing towards the end, like trapezoid self scheduling (TSS) andguided self scheduling (GSS). On heterogeneous COWs, GSS produced the most evenworkload. However, we observed that slowest processors became a bottleneck for the.
runtime, if they received the first (largest) batch of NTNP

50

5. Parallelized Quartet Puzzling

NR
Hence, we modified the GSS algorithm to schedule batch sizes of GSS(Ξ) = 2NP,
to reduce the possible waste of computation time. The modified GSS algorithm wasapplied to schedule the computational more expensive pairwise distance computations
in step (5.1a).

To apply GSS efficiently to the computationally most demanding parts of the program,
ML step and puzzling step, an additional adaptation was necessary. An atomic batch
size was introduced due to coding specificities of 4 bit per quartet to store the preferred
quartet topologies from the ML step (cf. 5.5), to make efficient transmission of theresults possible. The atomic size can also be used to keep the batch size at a certain
level to reduce communication when only a few tasks are left to schedule.
Based on the modified GSS, the adaptation used for distributing the tasks of the ML
step and the puzzling step determined the batch size by the formula
if NR < ABS ,

 NR


 

 NR 
 NR 

NR
+ABS
−
mod
ABS
if
mod
ABS
> 0 (5.1)
SGSS(Ξ) =
2NP
2NP
2NP




  NR 
otherwise.
2NP
The second line of Eq. 5.1 rounds to the next multiple of the atomic size ABS , ifnecessary. The first line schedules the rest NR if already smaller than the atomic size
ABS . The adaptation produces a smoothed distribution of batch sizes, therefore it isreferred to as smooth guided self scheduling (SGSS).
For ML step and puzzling step of the QP algorithm, we determined an atomic size of
4 tasks to be a good tradeoff between a desirably small atomic size and the reduction ofcommunication overhead.

5.7. Efficiency of Parallel TREE-PUZZLEWhen analyzing the efficiency of parallelized algorithms the scalability of the parallelized
code with increasing numbers of processors is relevant. The scalability was evaluated on
simulated and biological data for benchmarking.

5.7.1. Benchmark Datasets and Setup
To analyze the efficiency of the parallel implementation, we used all datasets describedin 5.3.
For each dataset parallel TREE-PUZZLE was applied to obtain benchmarks for 2,
4, 8, 16, 32, 48, and 64 processors on a Cray T3E-1200 parallel computer. Note, we
use processes and processors interchangeably because the queuing system guaranteed5.7. Efficiency of Parallel TREE-PUZZLE

51

that one processor executes exactly one process. Due to the local setup, the sequential
version of TREE-PUZZLE did not run with the same priority as parallel jobs. Thus,
the speedup is compared to a 2 processor run, i.e. 1 master and 1 worker process. The
abscissa in the diagrams of Figures 5.3 to 5.5 shows the number of worker processors.

5.7.2. Results for Parameter Estimation
Fig. 5.2 shows that the parallel implementation of the parameter estimation with themodified GSS algorithm in step (5.1a) and SC in step (5.1c) does not scale for large
numbers of processors. Although all calls of steps (5.1a) and (5.1c) accumulate over
90% of the running time (Tab. 5.1), every single call is very fast. The parallelized partis then interrupted by a sequential part. Due to this fine granularity, the result is not
surprising. However, a slight increase in scalability can be observed with respect to taxa
and sequence length. For the parallel estimation of parameters for datasets with n > 200
taxa the gain of runtime reduction is substantial.
The biological dataset with 50 taxa and 289 amino acids (aa) fits well with simulated
data for 48 taxa and 1000 nt. The high computation time of the shorter aa dataset is
attributed to the increased size of the alphabet (20 aa to 4 nucleotides).

5.7.3. ML Step and Puzzling Step
While the scalability of the parameter estimation is dissatisfactory, the parallel ML step
(Fig. 5.3) and the parallel puzzling step (Fig. 5.4) show an almost perfect speedup indicating a very good performance of the parallelization strategy. Concerning the biologicaldataset the speedup is almost perfect in the ML step, but decreases in the puzzling stepwith growing numbers of processors.
While the runtime for the biological dataset is highest in ML step, the sequencelength has no influence on the runtime of the puzzling step because only the quartet tree
topologies are used here.

5.7.4. Overall Scaling of Parallel TREE-PUZZLE
Despite the disappointing scalability of the parallel parameter estimation and the additional sequential consensus step, the overall runtimes of the TREE-PUZZLE runs stillshow convincing scalability (Fig. 5.5). The scaling increases with growing numbers oftaxa, and for 64 sequences of length 1000 nt the speedup is almost perfect. Although this
could not be shown with the protein dataset used increasing the number of sequences
should bring a better speedup for protein dataset as well.
As in the ML step and the parameter estimation, the computation time consumed by
the non-parallelized consensus step increases from DNA to amino acid data because of
the higher alphabet size (Fig. 5.3, Tab. 5.2).

52

5. Parallelized Quartet Puzzlinga.

Speedup in Parameter Estimation
70

50 seq, 289aa
32 seq, 500nt
32 seq, 1000nt
48 seq, 500nt
48 seq, 1000nt
64 seq, 500nt
64 seq, 1000nt
128 seq, 500nt
128 seq, 1000nt
perfect speedup

60

Speedup

50
40
30
20
10
0
0

10

20
30
40
50
Number of Worker Processesb.

60

Runtime in Parameter Estimation
10000

Runtime

1000

100

10

1
0

10

20

30

40

50

60

Number of Worker Processes

Figure 5.2.: Speedup (a) and runtime(b) of the parameter estimation for the different
datasets.

5.7. Efficiency of Parallel TREE-PUZZLEa.

53

Speedup in ML Step
70

50 seq, 289aa
32 seq, 500nt
32 seq, 1000nt
48 seq, 500nt
48 seq, 1000nt
64 seq, 500nt
64 seq, 1000nt
perfect speedup

60

Speedup

50
40
30
20
10
0
0

10

20
30
40
50
Number of Worker Processesb.

60

Runtime in ML Step
100000

Runtime

10000
1000
100
10
1
0

10

20

30

40

50

60

Number of Worker Processes

Figure 5.3.: Speedup (a) and runtime (b) of the maximum-likelihood step of Quartet
Puzzling for the different datasets.

54

5. Parallelized Quartet Puzzlinga.

Speedup in Puzzling Step
70

50 seq, 289aa
32 seq, 500nt
32 seq, 1000nt
48 seq, 500nt
48 seq, 1000nt
64 seq, 500nt
64 seq, 1000nt
perfect speedup

60

Speedup

50
40
30
20
10
0
0

10

20
30
40
50
Number of Worker Processesb.

60

Runtime in Puzzling Step
100000

Runtime

10000
1000
100
10
1
0

10

20

30

40

50

60

Number of Worker Processes

Figure 5.4.: Speedup (a) and runtime (b) of the puzzling step of Quartet Puzzling forthe different datasets.

5.7. Efficiency of Parallel TREE-PUZZLEa.

55

Overall Speedup
70

50 seq, 289aa
32 seq, 500nt
32 seq, 1000nt
48 seq, 500nt
48 seq, 1000nt
64 seq, 500nt
64 seq, 1000nt
perfect speedup

60

Speedup

50
40
30
20
10
0
0

10

20
30
40
50
Number of Worker Processesb.

60Overall Runtime
100000

Runtime

10000

1000

100

10
0

10

20

30

40

50

60

Number of Worker Processes

Figure 5.5.: Overall speedup (a) and runtime (b) of the parallelized TREE-PUZZLE
program for the different datasets.

56

5. Parallelized Quartet Puzzling

5.8. Discussion
To compute ML trees, it is necessary to estimate the underlying model parameters. This
chapter illustrates a first attempt to speed up the estimation procedure by parallelizing
the estimation heuristic implemented in TREE-PUZZLE. The combination of message
passing with SC and a modified GSS leads to reduced computation time. However, thespeedup is not even close to perfect. A factor of 17 for 32 processors will neverthelessenable biologists to study large datasets. For example, the estimation of evolutionary
parameters and rate heterogeneity parameters for 215 rRNA sequences of length 2100
nt took only 8.6 minutes on a homogeneous 32-processor COW instead of 2.5 hours ona single processor. This reduction in computation time makes TREE-PUZZLE also aworthwhile stand-alone application to estimate evolutionary parameters. Thus in the
future, sub-sampling strategies for parameter estimation (Meyer et al., 1999; Meyer andvon Haeseler, 2003) may no longer be necessary for moderately sized datasets. This isa first step to develop more efficient and scalable algorithms to support large phylogeny
reconstructions.
The amount of data available naturally leads to large tree reconstruction projects,
e.g., to reconstruct a tree for each of the 25,000 multi-sequence protein families in the
SYSTERS database (Version 3, Krause et al., 2002). For the smaller families the trees
can be reconstructed using the sequential TREE-PUZZLE. However, for large families(n ≥ 100 sequences) the parallelized version might be an alternative, because the presented parallelization provides for the first time a parallel tree reconstruction approach
for amino acid sequences.
In the past the DNAML program (Felsenstein, 1981) was subject to parallelization
attempts. Several parallel implementations were published showing moderate to wellscaling behavior (Olsen et al., 1994; Schmidt, 1996; Trelles et al., 1998; Stewart et al.,
2001; Stamatakis et al., 2002).
This chapter presented the parallel version of the QP algorithm which is an integral
part of the TREE-PUZZLE package (Schmidt et al., 2002; Strimmer and von Haeseler,
1996). The parallelized QP shows convincing scalability. One advantage of the parallel
TREE-PUZZLE is its simplicity. The master/worker setup with one master process
for the scheduling is sufficient to handle all tasks and communication. No additional
processes like foreman or dispatcher processes in the fastDNAml implementations (Olsenet al., 1994; Stewart et al., 2001) are necessary.
Since Quartet Puzzling is the only parallelized ML-based reconstruction method forprotein phylogenies so far, it is a major contribution for the ML analysis of large protein
datasets.
Although we have made considerable progress, more work needs to be done. The maindrawback in scalability remains with the merging of the split lists from the workers into
a consensus tree and the computation of ML branch lengths of the consensus tree. How
to parallelize the computation of ML branch lengths remains an open problem.

6. Large ML Trees from Sequences
Using Quartets
6.1. Introduction
As described before, the ML step and the puzzling step are the two steps demanding
most of the runtime in the Quartet Puzzling procedure (a detailed runtime analysis of
QP has been presented in 5.4). In chapter 4 we demonstrated that the complexity ofthe puzzling step can be reduced by one order of magnitude using an efficient O(n4 )algorithm, which leads to substantial speedup. However, the complexity Θ(n4 ) of theML step cannot be reduced, since the QP algorithm requires the evaluation of all 3 n4
possible quartet topologies. For datasets of 200 sequences or more the computation timerequired to evaluate all quartets is very large. For example, the evaluation of all quartet
topologies for a red algae dataset of 215 sequences from the European small subunit(ssu) rRNA database (Van de Peer et al., 2000b) assuming rate heterogeneity would
take more than 6,500 hours on a SUN Enterprise 450 server.
Hence, a strategy is suggested that does not require the computation of all quartet
trees dividing the set of sequences S into subsets. In order to ensure that our algorithm
works, a minimal amount of pair-wise overlap, i.e., the number of sequences shared by
the different subsets of the alignment is required. Once the quartet trees are computed,
we proceed with a modified version of the puzzling step (cf. 3.2.2). An overall consensus
tree is constructed from a series of so-called intermediate trees produced in the modified
puzzling step algorithm. We show, that the resolution of the consensus tree criticallydepends on the number of quartets actually evaluated.
The outline of the chapter is as follows, in section 6.2 we briefly introduce the relevant additional notation followed by the explanation of the modified Quartet Puzzling
algorithm for large datasets. The subsequent section gives some insights about the correlation between the resolution of the tree and the amount of computing time spend.
The final section discusses applications and gives an outlook on further modifications of
the proposed method.

6.2. The Modified Quartet Puzzling Algorithm
The Quartet Puzzling algorithm requires the computation of n4 quartets, which is often
not feasible if n is large. This chapter introduces a modified version of the algorithm

57

58

6. Large ML Trees from Sequences Using Quartetsby computing only quartets for a collection of subsets S1 , S2 , S3 , . . . , Sk of S. Then Qg
describes the set of all |S4g | quartets from Sg for all g = 1, 2, . . . k. We require thatk[

Si = S,(6.1)|Si | ≥ 4, for i = 1, . . . , k

(6.2)

i=1

( i
[

)
Sj

∩ Si+1≥ 3, for i = 1, . . . , k − 1.

(6.3)

j=1

We call equation 6.3 the overlap condition. It ensures that we can assemble an nspecies tree from the quartet trees. Note, that according to the overlap condition wewould have to compute

k 
X
|Si |

(6.4)

4

i=1

quartets only, dividing S into k subsets.

6.2.1. The Algorithm ModPUZZLE
Before discussing some computational aspects, we will explain the algorithm for k = 2.
If k > 2, then one has to apply the algorithm k − 1 times by setting

S1 =

j
[

Sj

(6.5)

S2 = Sj+1

(6.6)

i=1

for j = 1, . . . , k − 1. The order of the subsets is chosen randomly and that the overlap
condition is fulfilled.
The ModPUZZLE algorithm works as follows (cf. also Figure 6.1)
step 1: Assume that an intermediate tree T (S1 ) is reconstructed using the normal puzzling step algorithm (cf. 3.2.2 and 4.4).
step 2: Compute
S1∩2 = S1 ∩ S2 .
The elements of S1∩2 are already in T (S1 ).

6.2. The Modified Quartet Puzzling Algorithm

59

step 3: T (S1∩2 )step 1/2: T (S1 )root-nodes
sequences in S1 and S2 (S1∩2 )step 6:

step 5:

X

X



 



  

edges with minimal penalty (Emin )and adjacent rooted subtrees

new root-nodesdeleted root-nodes

Figure 6.1.: The modified Quartet Puzzling algorithm. For details see text.
step 3: Derive the subtree T (S1∩2 ) from T (S1 ), that has as leaf-set the sequences in S1∩2
and as internal node set all nodes from T (S1 ) that are on a path connecting two
leaves from S1∩2 . Each node of degree 2 in T (S1∩2 ) is then a root of a subtree from
T (S1 ) with leaves in S1 but not in S2 . Nodes of degree 2 are called root-nodes.
step 4: Pick a sequence X ∈ S2 \S1 and compute the penalties of each edge in T (S1∩2 )using the normal puzzling step algorithm (see 3.2.2 and 4.4).
step 5: Let E be the collection of all edges in T (S1∩2 ) with minimal penalty, then insert
X randomly either on one edge e ∈ Emin or on an edge from a rooted subtree ofT (S1 ) whose root node has at least two edges in T (S1∩2 ) with minimal penalty.
step 6: If X is the inserted sequence, then update T (S1∩2 ) by adding X to the leaf-set
and by adding all nodes of T (S1 ) that lie on the path between X and the root-node
R of the subtree of T (S1 ), where X was inserted, to the set of root nodes. Finally,
delete R from the root-node set.
step 7: Delete X from S2 \S1 . If S2 \S1 6= ∅ goto step 4, otherwise goto step 8.
step 8: Repeat step 1 to step 7 of the algorithm several times for different input ordersof the subsets and the sequences within the subsets.

60

6. Large ML Trees from Sequences Using Quartets

step 9: From this collection of trees the majority consensus tree is calculated.

6.3. Computational Aspects
To simplify matters, we assume that we partition the data in k ≤ n subsets of size ν = nk +o, where o defines the overlap between two neighboring sets, i.e. |Si ∩ S(i mod k)+1 | = o for
i = 1, . . . , k, where we require that o ≥ 3. Although not crucial for the ModPUZZLE
algorithm, for the practical application but to add freedom to the repetitions of steps
1 to 7, we now construct the last subset Sk overlapping S1 .
For k = n and o = 3 we have to compute n − 3 quartets only. This, however, is
not enough information to reconstruct the tree, even if the data were perfectly tree-like.
Thus, one should try to increase o as much as possible, guided by the time one is willingto spend for the calculation of thek  
ν−o
o
+
4
4

(6.7)

possible quartets, where o > 3.

6.4. Some Practical Measures on the Method
In phylogenetic tree reconstruction one wants to reconstruct a fully resolved tree, i.e., allinner nodes of the tree have degree 3. However, due to the lack of a clear phylogenetic
signal it is not always possible to achieve this goal (cf. 3.1). In these situations the tree
reconstructed is not fully resolved. The amount of resolution can be measured by the
number of splits, i.e. inner branches, that partition the sequences into two non-empty
subsets. To analyze the influence of the size of the overlap between subsets on the
resolution, we carried out a simulation study. To this end, we simulated the evolutionof DNA sequences on a tree with 50 leaf-vertices using the Seq-Gen package (Rambaut
and Grassly, 1997).
The ModPUZZLE algorithm was applied to the resulting data by randomly splitting
the 50 sequences into k = 5 subsets of varying size and different overlap. Table 6.1
summarizes the results. The results show that one should try to maximize the overlap,
since larger overlap enables the a higher resolution of the reconstructed tree. It showsthat an overlap of 20 induces the computation of roughly 50% of all possible quartets.
On the other hand, 39 from the 47 possible splits could be recovered, which provides areasonably good resolution.
The ModPUZZLE algorithm was also applied to the alignment of all 215 red algae ssu
rRNA sequences from the European small subunit rRNA database (Van de Peer et al.,2000b) as a biological dataset. Due to the large number of quartets (86,567,815 possible

6.5. Discussion and Possible Extensions

61

Table 6.1.: Results of the tree reconstruction using a simulated alignment of 50 DNA
sequences with increasing overlap.
kνo number of% of resolved % of
subsets sequences overlap
quartets quartets
splits splits
5
13
3
3,230
1.4
2
4.3
5
15
5
6,800
2.6
7
14.9
5
20
10
23,175
10.1
23 48.9530
20
112,800
49.0
39 83.0
full set50
230,300
100.0
47 100.0

quartets) tests were only ran for a limited set of values for k and ν. Nevertheless, thislarge biological dataset shows similar characteristics as the simulated one. The resolution
increases with the number of quartets available for the tree reconstruction.
The tests on the simulated data as well as the limited tests on biological data showthat with the minimal amount of overlap (o = 3) almost no resolution of the trees can begained. However, the resolution increases with the number of quartets used. Remarkable
is the effect that the percentage of resolved splits grows faster than the percentage ofquartets used (Tab. 6.1). Therefore it seems to be possible to reconstruct resolved trees
even if one does not use all quartet trees. The potential runtime reduction achieving
similar resolution, however, crucially depends on the amount of phylogenetic information
present in the alignment.

6.5. Discussion and Possible Extensions
In this chapter a very simple algorithm has been suggested to reconstruct phylogenetic
trees from large datasets. The method is based on a modified version of the puzzling
step algorithm described in 3.2.2 and 4.4. The modified algorithm has the flexibility toadjust the amount of computing time one is willing to spend. If one is only interested in
the coarse structure of the tree, then one needs to compute only very few quartets, thus
obtaining a more or less unresolved tree. If one wants the fine details of the ramificationsof the tree, one needs to compute a lot more quartets by increasing the overlap between
the subsets.
But, as shown in Table 6.1, the percentage of resolved splits seems to grow faster than
the amount of quartets used. This observation leads to a strategy how to analyze large
datasets, which is not fully exploited here. Instead of randomly assigning sequences to
the k subsets once in some kind of linear order, one could use more decompositions tobuild a network of subsets. To increase the resolution of the final tree, one may also use
a data guided approach. For example, a threshold graph (Barthelemy and Guenoche,
1991; Huson et al., 1999) based on the pairwise distances can be used as an indicator62

6. Large ML Trees from Sequences Using Quartets

how to group sequences. The applicability of this strategy needs to be analyzed bysimulations.
Another extension of the algorithm is also possible. Instead of analyzing the gene treeof one set of aligned sequences, we may very well assume that each subset S1 , S2 , . . . ,Sk contains the collection of species for which a sequence alignment is available for one
out of k different genes. We then can compute the tree for the entire set of species S,
based on k different genes, without requiring that one gene sequence is known for all n
species in S. The properties of this approach need to be investigated.
In chapter 7 a method will be described similar to the latter extension.

7. Phylogenetic Trees from Multiple
Genesets with Missing Data
7.1. Introduction
The large amounts of molecular sequence data currently available serve two needs inthe phylogenetic analysis of the relationship of species. On the one hand, the numberof interesting species available for analysis grows (vertical growth). On the other hand,more different sequences per species get available (horizontal growth). Unfortunately,
the number of sequences is not distributed evenly among species of interest. This often
makes collecting a dataset for phylogenetic analysis a painful decision on the trade-off
between the amount of taxa and the number of sequences.
While methods abound to infer phylogenies from a set of aligned sequences (Swofford
et al., 1996), only a small number of methods exists to combine data of different genes,
proteins, or genomic areas for joint analysis.
There has been a large ongoing debate how to combine different datasets to reconstruct
trees (cf. de Queiroz et al., 1995, Bininda-Emonds, 2002, and Page and Holmes, 1998,
chap. 8 for review). Two paradigms have been discussed, which we will classify by the’distance’ of the combination event from the underlying data into low level and highlevel methods (cf. Fig. 7.1). Note that low and high does not imply a quality rating.
The first paradigm is total evidence (also often called combined or simultaneous analysis) where the data is combined directly by concatenating the alignments. Hence, total
evidence methods will be referred to as low level methods (cf. Fig. 7.1).
The other paradigm is the so-called separate analysis (also called taxonomic congruence and consensus or supertree approach). Such methods combine trees reconstructed
separately from the single datasets. Since the combination takes place far from theunderlying data, such methods are referred to as high level methods (cf. Fig. 7.1).
Both paradigms have their advantages and drawbacks (for review see de Queiroz et al.,
1995, Bininda-Emonds, 2002, and references therein). Major criticisms are, e.g., theproblems to jointly model the evolutionary process for the concatenated dataset in low
level methods and the loss of information in the high level methods since the underlying
data is in general not considered when the trees are combined. These arguments gave
rise to the design of an alternative method which will be proposed in this chapter.
The method presented here follows a medium level paradigm, according to the level ofcombination (cf. Fig. 7.1). Before particularizing the procedure (section 7.3), commonly

63

64

7. Phylogenetic Trees from Multiple Genesets with Missing Data

used high and low level methods are described (section 7.2). Then the new method isapplied to a biological dataset in section 7.4. The results are discussed in the context
of total evidence and consensus/supertree methods. Finally, problems will be discussed
and an outlook on further improvement and extensions will be given.

7.2. Methods to Combine Datasets
7.2.1. Low Level Combination: Total Evidence
As mentioned above, low level methods are commonly referred to as total evidence,
combined or simultaneous analysis methods. Claiming that all information (evidence)
available should be used for phylogenetic analysis, all single datasets are combined into
one single ’supermatix’ (Sanderson et al., 1998) as shown in Fig. 7.1. This is achievedby concatenating all source alignments filling missing data with gap characters. The
overall alignment is then used as input for phylogenetic analysis. This method seems
unproblematic for complete datasets, where for each species one sequence in each source
dataset exists. How total evidence approaches handle missing data crucially depends onthe method used to reconstruct a phylogeny from the concatenated alignment. Sincesome programs exclude alignment columns with gaps, such programs will end up onlywith data from genes which are available for the whole set of species. Datasets withmissing sequences are discarded.

7.2.2. High Level Methods
While low level methods combine the source datasets directly, high level methods construct a tree for each source dataset. The trees from the source datasets are then combined to an overall tree (Fig. 7.1). If all source datasets contain the full set of species,
consensus techniques can be applied combining all the equally sized trees into one consensus tree. To combine sets of trees with unequal, but overlapping sets of species,
so-called supertree methods have been developed which amalgamate the input trees into
one overall supertree.
7.2.2.1. Combining Equal-Sized Datasets: Consensus Methods
Consensus methods aim to construct a consensus tree from a set of trees preserving
common topological details of the source trees. Commonly applied consensus methods
are strict consensus, majority-rule consensus (Margush and McMorris, 1981), semi-strict
consensus (Bremer, 1990), and Adams consensus (Adams III, 1986).
Most consensus methods are based on splits or bipartitions (cf. 2.2.1) found in a
majority of the source trees (Ml consensus sensu McMorris and Neumann, 1983, where

7.2. Methods to Combine Datasets

65

low level

combination

























method










































































































low levelhigh levelhigh levelhigh levelhigh levelmedium level

Total Evidence













medium level




































high level







































































splitsConsensus











































































































MRP, MRF







































ModMinCut,
MinCutSupertree












































































Quartet-based
supertree







































































Medium levelcombination

























Figure 7.1.: Level of combination with regard to distance from the underlying datasets
l represents the percent occurrence of included splits). Sets of non-contradicting bipartitions are chosen to reconstruct the consensus tree.
The majority-rule consensus (Margush and McMorris, 1981) uses all splits occurring in
the majority of input trees. Usually a majority-rule consensus is assumed (Ml consensus
with l > 0.5, McMorris and Neumann, 1983). This means all splits occurring in morethan 50% of the source trees (see Fig. 7.2). Considering l > 0.5 ensures that all splits

66

7. Phylogenetic Trees from Multiple Genesets with Missing DataA

D
AB|CDEF

B

E

A
strict consensus
ABC|DEF
B

D

C

F

C

D

A

E

B

F

A

D

C

E

B

F

E

ABC|DEFCF

Tree AAD

ABC|DEF − 3 (100%)
AC|BDEF − 2 (66.7%)
AB|CDEF − 1 (33.3%)
ABCD|EF − 1 (33.3%)

AC|BDEFC

E
ABC|DEF

B

F

Tree B

A

D

C

E

B

AC|BDEF
ABC|DEF
ABCD|EFF

Tree Csemi−strict
ABC|DEF
ABCD|EFmajority−rule
AC|BDEF
ABC|DEF

Figure 7.2.: Consensus methods
A

B

C

Tree A

D

E

E

B

C

Tree B

D

A

A

B

C

D

E

strict/majority consensus

B

C

D

A

E

Adams consensus

Figure 7.3.: Adams Consensus
can be incorporated in a tree topology, because no two contradicting splits are able tooccur in more than 50% of the input trees.
For the strict consensus tree only bipartitions are chosen that are found in all input
trees (M1.0 , see Fig. 7.2). The semi-strict consensus (Bremer, 1990) contains all splitsthat are not contradicted in any of the input trees, e.g. the split ABCD|EF is uncontradicted by the ABC-trifurcation in tree A and B in Fig. 7.2. Note that such splits
can occur in less than half of the source trees. If all trees are binary trees strict andsemi-strict consensus will always produce the same trees.
Contrary to most consensus methods which are based on bipartitions and their percent
occurrence in the source trees, the Adams consensus is based on common nestings intrees, i.e., taxa frequently occurring together in the same subtree. This method tries
to find groups of sequences that are commonly occurring together in the source trees(Adams III, 1986). Adams consensus trees can contain groups that cannot be found inany of the source trees. This makes it difficult to interpret. Yet it can be informative
when there are sequences that are difficult to place. Such sequences are moved to theroot of the Adams consensus tree (Fig. 7.3).
7.2.2.2. Combining Overlapping Datasets: Supertree Methods
Supertree methods have been developed to combine sets of overlapping trees into ’supertrees’ containing all leaves found in the source trees. Some of the methods are re-7.2. Methods to Combine Datasets
A

B

C

1

D

67
E

F

A

E

G

3

H
6

5
2
4
Ragan/Baum

A
B
C
D
E
F
G
H

Purvis

1 2 3 4

5 6

1
1
0
0
0
0
?
?

0
?
?
?
1
?
1
1

1
1
1
1
0
0
?
?

0
0
1
1
0
0
?
?

1
1
1
1
1
0
?
?

0
?
?
?
0
?
1
1

A
B
C
D
E
F
G
H

1 2 3 4

5 6

1
1
0
0
?
?
?
?

0
?
?
?
1
?
1
1

1
1
1
1
0
?
?
?

0
0
1
1
?
?
?
?

1
1
1
1
1
0
?
?

?
?
?
?
0
?
1
1

Figure 7.4.: Coding schemes to encode tree topologies for MRP methods: Baum/Ragan
and Purvis’ scheme
stricted to rooted trees, based on the argument that the broad majority of trees publishedare rooted trees. Input trees cannot be produced from different datasets only, but canalso be obtained from the literature.
Matrix Representation Methods A commonly used method to construct supertrees
is Matrix Representation using Parsimony (MRP, Baum, 1992; Ragan, 1992), whichuses binary coding to represent the input trees. The splits induced by all source trees
are coded into a matrix with binary characters (’0’, ’1’) with missing data (’ ?’). The
binary matrix is then used to construct an overall tree (cf. Fig. 7.1) using Maximum
Parsimony methods (Swofford et al., 1996). The most abundant coding schemes arethose independently suggested by Baum (1992) and Ragan (1992) and the modified
scheme by Purvis (1995). Baum/Ragan code the bipartitions of an input tree assigning’1’ to all leaves in the one part of the bipartition and ’0’ to the other that contains theroot. All missing taxa are assigned the missing data character ’ ?’ (Fig. 7.4). Purvis(1995) differs in that only the sistergroup of a clade is assigned ’0’, while all other leavesof that part outside the clade are assigned ’ ?’ (Fig. 7.4).
To give input trees different weights in the analysis, different weighting schemes havebeen suggested (see Salamin et al., 2002, Sanderson et al., 1998, and references herein
for more details).
Recently, another method has been proposed called Matrix Representation using Flipping (MRF, Chen et al., 2003), which also uses a binary matrix coding of the trees
as described above. The MRF supertree is constructed by ’flipping’ conflicting character states from ’1’ to ’0’, or vice versa to produce a matrix that does not containcontradictions. Their optimal solution is the tree that needs the least ’flips’.

68

7. Phylogenetic Trees from Multiple Genesets with Missing Data

Direct Supertree Methods Another family of supertree methods use a more graphtheoretical approach. Members of that family are the OneTree/BUILD algorithm(Aho et al., 1981; Bryant and Steel, 1995), MinCut Supertree algorithm (Semple
and Steel, 2000), and the modified MinCut Supertree algorithm (ModMinCut
Supertree, Page, 2002). These algorithms take as input sets of rooted trees. The
BUILD algorithm and the OneTree algorithm produce a supertree only if the input
trees are compatible, i.e., the trees are not contradictory (cf. strict consensus in 7.2.2.1).
Since this is almost never the case for biological data, these algorithms are mainly usefulto test for compatibility of trees.
MinCut Supertree tries to yield a supertree applying a minimum cut algorithm tothe graph ST constructed from the set of input trees T , following Listing 7.1 (see also
Fig. 7.5 and Fig. 7.1). The notations follow Page (2002) – for further details see there.
Page (2002) suggested a modification (here called ModMinCut) to choose the minimum cuts, that takes into account preserving uncontradicted information from thesource trees. The ModMinCut Supertree algorithm marks all edges in the graph
ST (step 7.1a) not contradicted by any tree in T . If ST has to be cut in step 7.1b
minimum cuts are preferred which do not cut uncontradicted edges. Hence, the amountof uncontradicted information from the source trees is preserved in the supertree (cf.
Fig. 7.6).
Quartet-Based Supertree Recently a quartet-based supertree method has been proposed by Robinson-Rechavi and Graur (2001). Instead of using a matrix representation,
all source trees are decomposed into sets of quartets contained in the topologies (Fig. 7.1).
Each quartet is weighted by the highest support value found on the path of the quartet’s
middle edge in any of the input trees (the authors used the TREE-PUZZLE package tocompute the trees and support values). From this set of weighted quartets an overall
tree is constructed with the quartet method implemented in the AllTree program (Bendor et al., 1998). In a first step, this method computes all subsets of taxa and keeps
those satisfying the most quartets. In a second step, it performs an exact search on eachsubset constructing an overall tree satisfying most quartets.

7.3. Medium Level Combined Phylogenetic Analysis
7.3.1. Notation
In the present chapter we assume a collection S = {s1 , . . . , sn } of n species and k different
genes to be used for combined analysis.
With S1 , S2 , . . . , Sk we denote subsets of S, such that
|Sg | ≤ 4

for each

g

and

k
[
g=1

Sg = S.

(7.1)

7.3. Medium Level Combined Phylogenetic Analysis
A

B

C

D

Tree A

E

A

B

C

69

D

Tree B
A

B

C

D

E

B
2

A,B

1
C

A

1
1

C

1

1

Supertree

1
D

E
S {Tree A,B}

1

D

E
Max

S {Tree A,B} /E {Tree A,B}

Figure 7.5.: The MinCut Supertree Algorithm. For more details, see text.
Listing 7.1: MinCut Supertree (Set of rooted trees T , set of n leaves)switch number of leaves n docase (n = 1) return a single node x1 ;case (n = 2) return a tree with 2 leaves x1 and x2 ;otherwise
7.1aConstruct the graph ST by connecting each pair of taxa a and b by an edge
e(a,b) if they are connected in at least one tree in T by a path that does notpass the root of that tree, i.e., MRCA(a, b) 6= root (cf. Fig. 7.5, S{TreeA,B} );Weight the edge e(a,b) by the number of trees supporting MRCA(a, b) 6= root;

if ST is a disconnected graph then
Let Si be the components of ST ;else
Construct ST /ETM ax as follows: Merge all nodes a and b connected by an edge
e(a,b) with maximum weight |T |, i.e., supported by the whole set of trees andremove edge e(a,b) (cf. Fig. 7.5);
7.1bRemove all edges from ST /ETM ax that are a minimum cut set of ST .
Let Si be the resulting components of ST /ETM ax ;foreach component Si do
Ti = MinCut Supertree (T |Si , L = Si ), where T |Si denotes the trees fromT with all leaves r 6∈ Si pruned;
Construct a new tree T by connecting the roots of Ti to a new tree;
return T

70

7. Phylogenetic Trees from Multiple Genesets with Missing Data

C
B
A
Y1
Y2
Y3
Y4

A
B
C
X1
X2
X3

Tree 1

Tree 2

C
X1
X2
X3
B
A
Y1
Y2
Y3
Y4
MinCut Supertree

A
B
C
X1
X2
X3
Y1
Y2
Y3
Y4
ModMinCut Supertree

Figure 7.6.: The ModMinCut Supertree algorithm. For more details, see text.
Each Sg represents the subset of species for which a multiple alignment based on geneg is available. These subsets will be called genesets. Finally, Qg , g = 1, . . . , k represent
the corresponding sets of possible quartets.
Instead of reconstructing the trees T (Sg ) for each set Sg , we will compute an overall
tree T (S) combining the information provided from the evaluation of the quartet sets
Q1 , Q2 , . . . , Qk from the k different alignments. Note that a taxon is not represented by
one sequence any more.

7.3.2. The Combined Quartet Method to Combine Genesets
7.3.2.1. Combining the ML Quartets
In the method proposed, the genesets will be combined on the level of the quartets. As
(g)
(g)
(g)
a guideline to combine the quartets we use the log-likelihood `ab|cd , `ac|bd , and `ad|bc forquartet {a, b, c, d} ∈ Qg on the geneset Sg .
To combine the datasets for each geneset S1 , S2 , . . . Sk , all three log-likelihoods foreach quartet in Q1 , Q2 , . . . Qk are evaluated first (cf. ML step in 3.2.2).
Then, we compute`ab|cdkX(g)=`ab|cd(7.2)

g=1`ac|bdkX(g)=`ac|bd
g=1(7.3)

7.3. Medium Level Combined Phylogenetic Analysis`ad|bc =

71

k
X
(g)
`ad|bc

(7.4)

g=1
(g)

for each quartet. For the sake of simplicity, the log-likelihoods `τwith the sequences L(τ ) is not represented by alignment g.

= 0 if the quartet

Alternatively, one can also compute`ab|cd| {g : {a, b, c, d} ∈ Qg } |`ac|bd=| {g : {a, b, c, d} ∈ Qg } |`ad|bc.
=| {g : {a, b, c, d} ∈ Qg } |`ab|cd =(7.5)`ac|bd(7.6)`ad|bc(7.7)`τ can be viewed as the average support a quartet tree τ receives from the set of sequencealignments. In case that a quartet is not represented by any alignment the averages are
set equal to zero.
7.3.2.2. Computing the Overall Tree
To reconstruct a phylogeny T (S) based on the collection of log-likelihoods from Equations 7.2 to 7.4 or 7.5 to 7.7 we apply the PUZZLE idea from 3.2.2. The straightforward
application of the QP algorithm is only possible if information exist for each quartet in
Q. This, for instance, is the case if at least one alignment Sg comprises all the species
in S.
If some quartets are missing, which usually happens if each Sg is a proper subset ofS, then these quartets are treated as unresolved. A quartet tree is selected randomly
among the three possible topologies. In this case it is necessary to examine whether theoverlapping genesets can be combined. This will be explained in the following.
7.3.2.3. Assessing Whether Genesets Can Be Combined
To combine two genesets Si and Sj a minimum pairwise overlap of 3 sequence among
the two subsets is required,|Si ∩ Sj | ≥ 3.

(7.8)

We call Eq. 7.8 the pair-overlap condition. The pair-overlap condition is different fromthe overlap condition (Eq. 6.3) in section 6.2 where the subsets were constructed fromone complete set of sequences to ensure combinability of the subsets. Here the subsets

72

7. Phylogenetic Trees from Multiple Genesets with Missing Data

are genesets. Their sizes and overlap are defined by the availability of sequence data forthe species in S. Note, that not every pair of genesets might share a sufficient overlap.
To assess whether these gene datasets can be combined, we construct an overlap graph
Govl = (V, E) with a set of nodes V (genesets) and a set of undirected edges E (overlaps)
with
V = {S1 , S2 , . . . Sk }E = {e(Si ,Sj ) for |{Si ∩ Sj }| ≥ 3} 1 ≤ i < j ≤ k.

(7.9)
(7.10)The edge weights consist of the amount of overlap fulfilling the pairwise-overlap condition(see for an example Fig. 7.7). This representation of the datasets provides insight into
several properties of the entire dataset.
If the graph consists of unconnected subgraphs, not all data is suited for a combined analysis. According to the pairwise-overlap condition gene datasets from different
connected components cannot be used simultaneously in the same analysis because noquartet information is are available to reasonably guide the insertion of the sequences.
Genesets from one connected component, however, can be used for combined analysis. Although genesets might exist within a connected component which do not sharesufficient overlap with each other to be combined, they can be linked by means of othergenesets which first have to be added satisfying the pairwise overlap condition, to connect
the two.
Note, that by applying the overlap condition from Eq. 6.3, sets from one connectedcomponent can gain an overlap ≥ 3 to the leafset L(Ti ) of a tree Ti so far reconstructedfrom genesets of another connected component. Nevertheless, combining those sets willproduce a very doubtful tree, because too little information is available to guide theinsertion of sequences.
The combinability of the genesets from a connected component is characterized by two
features of the overlap graph. A high connectivity within a component (each geneset
shows overlap to many other geneset) results in network-like quartet information among
the genesets. This reduces the need for mediating sets to connect non-overlapping genesets. Even more important is the size of the overlap. The higher the overlap, the morequartets are shared among two neighboring genesets. Consequently, a higher amount of
information is available to guide the insertion of sequences from a geneset into a tree
already reconstructed.
7.3.2.4. Overlap-Guided Puzzling Step
As described above, sufficient overlapping information is needed to reasonably insert a
new leaf si+1 into a tree Ti . Hence, instead of using a random permutation of leaves, thepermutation has to follow certain restrictions. To satisfy the pairwise overlap condition,
we apply a graph based procedure similar to Prim’s minimum spanning tree (MST)7.4. The Phylogeny of the Grasses

73

algorithm (see, e.g., Cormen et al., 2001) assuming equal edge weights. We start withthe leaves from a randomly picked geneset S 0 ∈ {S1 . . . Sk }. We define a front set F
containing all unused nodes (genesets) from the overlap graph that are connected by anedge to any geneset already used in the tree. From the sequences in S 0 we construct
the tree T (S 0 ) applying the usual puzzling step algorithm (see section 3.2.2 and chapter
4). Then we randomly draw one geneset S 00 from F. We remove S 00 from F and add
the sets connected to S 00 , but not yet used to F. The sequences of S 00 are then added
to the tree. Guided by the overlap graph we thus construct an intermediate tree by
sequentially adding the genesets. As in the usual puzzling step many intermediate trees
are constructed using different orders of leaves and genesets.
7.3.2.5. Relative Majority Consensus
Since missing data adds more ambiguity, this naturally hampers the construction of aresolved majority-rule consensus tree. We therefore decided to apply a majority consensus which is slightly different from the majority-rule consensus used in section 3.2.2.
Instead of using only splits occurring in more than 50% (Ml with l > 0.5) of the intermediate trees, also splits below 50% are considered. The aim is to use as many splits as
are supported by relative majority to extract a maximum uncontradicted informationfrom the intermediate trees.
The relative majority consensus Mrel adds all congruent splits from the set of intermediate trees in descending order of occurrence. No splits will be accepted for theconsensus that have the same or lower percentage occurrence as any incongruent split.
This consensus procedure does not imply a fixed threshold but uses a variable percentage
down to the first incongruence guided by the relative majority.

7.4. The Phylogeny of the Grasses
In the following, the grass (Poaceae) dataset provided by the GPWG (Grass Phylogeny
Working Group, 2001, http://www.ftg.fiu.edu/grass/gpwg/) is reanalyzed with the
new medium level method as well as with other methods for combined analysis.

7.4.1. The Dataset
We used all molecular sequence data from this dataset. There were three nuclear loci,
NADH dehydrogenase, subunit F (ndhF ), the internal transcribed spacer (ITS) of ribosomal DNA, and granule bound starch synthase I (GBSSI or waxy), as well as three chloroplast genes, ribulose 1,5-bisphosphate carboxylase/oxygenase, large subunit (rbcL), RNApolymerase II, β” subunit (rpoC2 ), and phytochrome B (phyB ). The morphological andrestriction data was omitted. The GPWG dataset consists of ’normal’ taxa as well asof composite taxa represented by sequences from several species of one genus (30 taxa)

74

7. Phylogenetic Trees from Multiple Genesets with Missing Dataor several genera (6 taxa). The overall number of taxa in the analysis is |S| = 66. The
available sequences are listed in Tab. 7.1.

7.4.2. Methods
The new medium level approach is compared to other high level as well as low level
methods. To this end, several tree reconstructions were performed using.
7.4.2.1. Low Level Combination
For total evidence analysis the alignments from the six genesets were concatenated into
one large supermatrix. From this large alignment an overall tree was constructed with
the TREE-PUZZLE program assuming the HKY model (Hasegawa et al., 1985). The
missing model parameters were inferred from the large alignment.
7.4.2.2. High Level Combination
For high level methods, first ML trees were reconstructed using the fastDNAml program(Olsen et al., 1994; Stewart et al., 2001). From the resulting set of six phylogenetic trees
supertrees were constructed applying a number of high level methods.
Two MRP supertrees were built using the Baum/Ragan encoding as well as Purvis’
scheme. The input trees were encoded using the SuperTree program (Version 0.8b)
by Nicolas Salamin. From the two matrix representations the most parsimonious trees
were constructed using the PARS program (see below). Equally parsimonious trees werecombined to a strict consensus tree (M1.0 ) with the CONSENSE program (see below).
In addition, supertrees were constructed with the MinCut (Semple and Steel, 2000)
as well as the ModMinCut Supertree algorithm (Page, 2002). The reconstruction
was performed by using Rod Page’s SUPERTREE software (Version 0.2) via the webinterface available at http://darwin.zoology.gla.ac.uk.
The choice of methods was guided by the availability of implementations or web interfaces. The programs CONSENSE and PARS belong to the PHYLIP package (Version
3.6a3, Felsenstein, 1993).
7.4.2.3. Medium Level Combination
The combined analysis was performed as described in 7.3.2. For each of the six genesets
all possible quartet trees were evaluated. The HKY model (Hasegawa et al., 1985) was
assumed. Branch lengths were optimized in the computation of the ML quartets (’exact
quartet’ option in TREE-PUZZLE). Missing model parameters were inferred from thecorresponding genesets. The resulting log-likelihoods were combined. The combined
likelihoods were then used to determine the set of supported topologies. Finally, the set

7.4. The Phylogeny of the Grasses

75

s
g

X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X

Avena gBromus s
Triticum gAristida s
Stipagrostis
Amphipogon s
ArundoMolinia g
Phragmites
Merxmuellera m.
Karoochloa
Danthonia sAustrodanthonia
Merxmuellera r.
Centropodia
Eragrostis s
UniolaZoysia s
Distichlis
Pappophorum g
Spartina s
Sporobolus s
Eriachne s
Micraira s
Thysanolaena
Gynerium
Chasmanthium s
Zeugites
Danthoniopsis s
Panicum s
Pennisetum s
Miscanthus g
Zea

X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X

X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X

X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X

X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X

- X
- X
X X
- X
- X
- X
- X
- X
- X
X X
X X
- X
X X
X X
X X
- X
- - - - - X
- X
- X
- X
- X
- X
- X
- X - X
X X
X X
X X

composite taxon, represented by sequences from several species
composite taxon, represented by sequences from several genera

ITS

X
X
X
X
X
X
X
X
X
-

waxy

ITS

X
X
X
X
X
X
X
X
X
-

rpoC

waxy

X
X
X
X
X
X
X
X
X
X
X
X
X
-

rbcL

rpoC

X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X

phyB

rbcL

X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X

ndhF

phyB

Flagellaria
Elegia s
Baloskion
Joinvillea s
Anomochloa
Streptochaeta s
Pharus s
GuaduellaPueliaEremitis
Pariana
Lithachne sOlyra sBuergersiochloa
Pseudosasa g
Chusquea s
Streptogyna
Ehrharta s
Oryza
Leersia s
Phaenosperma
Brachyelytrum
LygeumNardusAnisopogon
Ampelodesmos
Stipa s
Nassella s
Piptatherum s
Brachypodium s
Melica s
Glyceria s
Diarrhena s

ndhF

Table 7.1.: Sequences available for analysis per taxon in the Poaceae dataset (Grass
Phylogeny Working Group, 2001)

76

7. Phylogenetic Trees from Multiple Genesets with Missing Data

Table 7.2.: The Poaceae dataset; estimated parameters and quartet statistics
sequence origin# sequencesalignment length
constant sites 1
A content
C content
G content
T content
test failed2
ts:tv3 parameterts:tv3 S.E.
ts:tv3 ratio
Y:R ts ratio4
average distance# quartets
resolved
partlyunresolved
missing

ndhF
chloroplast
65
2210
10.1%
27.3%
16.5%
17.6%
38.7%
0
2.14
0.08
1.93
1.33
0.077
67704094.41%
5.00%
0.60%

phyB
nucleus
40
1182
5.6%
21.5%
26.9%
29.2%
22.3%
4
1.87
0.08
1.85
0.96
0.160
91390
92.91%
5.77%
1.33%

rbcL
chloroplast
37
1344
45.3%
27.1%
19.3%
24.9%
28.7%
0
1.69
0.12
1.66
0.82
0.053
66045
94.57%
4.98%
0.45%

rpoC
chloroplast
34
777
2.4%
40.5%
14.6%
28.2%
16.7%
1
2.31
0.23
2.98
0.21
0.107
46376
71.49%
9.00%
19.52%

waxy
nucleus
19
773
54.2%
21.3%
29.4%
33.8%
15.5%
2
1.11
0.08
1.05
0.63
0.147
3876
94.09%
5.62%
0.28%

ITS
nucleus
47
322
48.4%
18.7%
31.6%
33.1%
16.6%
2
1.58
0.11
1.45
0.85
0.154
17836585.98%
10.23%
3.79%

combined
–
–
–
–
–
–
–
–
–
–
–
–
–
–
72072091.87%
4.23%
0.08%
3.82%

1 including gapped positions

3 transition:transversion

2 sequences failing χ2 test on composition

4 pyrimidine:purine transition ratio

tot. evid.
mixed
66
6608
0%
26.4%
20.0%
23.0%
30.6%
43
1.77
0.04
1.73
1.01
0.087
72072094.77%
4.17%
1.06%

of quartet topologies and the overlap graph were fed to the modified TREE-PUZZLE
program to reconstruct the combined tree. 100,000 puzzling steps were performed tobuild an Mrel consensus (cf. section 7.3.2.5).

7.4.3. Parameter Estimates from Poaceae Dataset
The Poaceae dataset used for combined phylogenetic reconstruction consisting of ndhF,
phyB, rbcL, rpoC, waxy, and ITS has a total of 6608 aligned positions. The different
genesets show substantial differences in their features as well as in the estimated parameters (cf. Fig. 7.2). The number of species in the datasets ranges from 19 in waxy to 65
in ndhF (cf. also Fig. 7.7). The sequences are closely related resulting in a low average
distance from 0.053 (rbcL) to 0.160 (phyB ) and a high percentage of constant sites in thecase of ITS (48.4%), rbcL (45.3%), and waxy (54.2%). The nucleotide frequencies differ
strongly between the genes. The estimated parameters for the HKY model (Hasegawa
et al., 1985) also vary a lot. The transition:transversion ratios range from 1.05 (waxy)
to 2.98 (rpoC ), while the transition ratio between pyrimidines and purines varies from
0.21 (rpoC ) to 1.33 (ndhF ). Furthermore, the fact only 23 (35%) of the 66 sequencesin the concatenated dataset pass the χ2 test on the sequence composition (Fig. 7.2)
shows that the overall character frequencies might be inappropriate for a majority ofthe sequences. This as well as the variation of the parameters clearly demonstrates thatfinding one single parameter set that matches the joint dataset can be problematic when
doing ’simultaneous analysis’.

7.4. The Phylogeny of the Grasses

77

rbcL
37
36

26
29

12

rpoC

21

waxy

15

19

40

ITS
47

21

14

30

phyB
40

ndhF
rbcL
rpoC
ITS
phyB
waxy

•

36 34 46
• 26 29
• 31
•

waxy

31

ITS

19

rpoC

10

phyB

34

34
46

rbcL

65

ndhF

ndhF

40
21
21
30
•

19
12
10
15
14
•

Figure 7.7.: Overlap graph of the gene datasets: Line width is chosen to represent the
number of sequences in the overlap or the dataset, respectively. The tableshows the amount of pairwise overlap.
The quartet resolution of the genesets, i.e., the amount of fully resolved quartets, ishigh (up to 94.57%), except for rpoC, where only 71.49% of the possible quartets could
be resolved completely.
The resolution of the combined quartet set with 91.87% completely resolved quartets
was rather high. This shows that the combination of the datasets did not lead to a
substantial reduction of the quartet resolution. Yet a high number of resolved quartets
does not exclude contradiction among the inferred quartet topologies. Missing data
among the overlapping genesets induces a total of 3.82% of the 720,720 quartets forwhich no phylogenetic information is available.
The parameters estimated from the concatenated alignment (Tab. 7.2) also present
substantial differences compared to the parameters inferred from the separate genesets.
While the average distance among the concatenated sites is still low (0.087), the constantsites have vanished. This is attributed to the gaps introduced to the alignment to fill up
the genesets to 66 sequences.

7.4.4. Combinability of the Poaceae Dataset
The overlap graph in Fig. 7.7 was constructed from the overlap among the genesets which
are listed in Tab. 7.1. The line widths represent the amount of sequences in the genesets(frames) or the amount of overlap between two adjacent genesets (edges), respectively.
The overlap graph consists of one connected component. This means that all genesets
can be used together in a combined analysis. The connectivity of the graph is high, every

78

7. Phylogenetic Trees from Multiple Genesets with Missing Data

geneset is connected to every other geneset by an overlap > 3. In general, the overlapof any two genesets is higher than half the size of the smaller geneset. This means that
more than 50% of the taxa contained in a geneset is also contained in a bigger geneset.
According to the pairwise-overlap condition and the high connectivity of the graphthe genesets can be combined. Note, that this judgment only takes into account thephysical overlap of the genesets and, thus, determines only the technical combinabilityof the dataset.

7.4.5. Reconstructed Poaceae Phylogeny
The taxonomy and the nomenclature follow the Grass Phylogeny Working Group (2001).
For more details about the phylogeny of the grasses refer to the article.
7.4.5.1. Total Evidence Tree
The total evidence tree of the grasses (Poaceae) based on the concatenated GPWG
dataset using TREE-PUZZLE for the simultaneous analysis is displayed in Fig. 7.8. The
tree is rooted by the outgroup taxa Flagellaria, Baloskion, Elegia, and Joinvillea. Thenthe so-called early branching taxa branch off the tree. First the Anomochloid subfamily
which is shown to be a monophyletic sister group of Pharus. Next, the monophyletic
Puelioid subfamily forms a sister group to the main group of the grasses. Within thegrasses six subfamilies establish the PACCAD clade. The PACCAD subfamilies Aristidoids, Danthonioids, and a main part of the Panicoids could be recovered. Aristidoids
and Danthonioids together form a monophyletic group.
The three other grass subfamilies, Bambusoids, Erhartoids, and Pooids, form monophyletic subtrees. These three subfamilies are reported to form a monophyletic group,
the so-called BEP clade (Grass Phylogeny Working Group, 2001) for different datasets.
The BEP clade itself, however, could not be revealed.
Two other total evidence grass trees based on the MP criterion have been published(Grass Phylogeny Working Group, 2001; Salamin et al., 2002) and will be presented
later in section 7.4.6.
7.4.5.2. MRP Supertrees
The two MRP supertrees present different topologies (Figs. 7.9 and 7.10). The firsttree was reconstructed from a matrix encoded according to Baum/Ragan. 29 mostparsimonious tree topologies have been found and joined by strict consensus (cf. Baum,
1992; Ragan, 1992). The resulting tree (MRP-BR) could establish most of the grass
subfamilies (Fig. 7.9). Starting from the root, the early branching taxa branch off in
the order Anomochloids, Pharus, and then Puelioids. Within the grass subtree, the7.4. The Phylogeny of the Grasses

79

Flagellaria indica

Anomochloa marantoideai / Anomochlooid

65

Streptochaeta / Anomochlooid

56
57
95

Pharus / Pharoid
Guaduella marantifolia / Puelioid

53

Puelia ciliata / Puelioid
Aristida / Aristidoid/PACCAD

94

Stipagrostis zeyheri / Aristidoid/PACCAD

60

86
43
54
63

59

Phragmites australis / Arundinoid/PACCAD
Arundo donax / Arundinoid/PACCAD
Amphipogon / Arundinoid/PACCAD
Thysanolaena maxima / Centothecoid/PACCAD
Chasmanthium / Centothecoid/PACCAD
Zeugites pittieri / Centothecoid/PACCAD

69

Danthoniopsis / Panicoid/PACCAD
Panicum / Panicoid/PACCAD
Pennisetum / Panicoid/PACCAD

79
48

Miscanthus g. / Panicoid/PACCAD

74

Zea mays / Panicoid/PACCAD
Merxmuellera rangei / Chloridoid/PACCAD

88

70

40

Centropodia glauca / Chloridoid/PACCAD
Eragrostis / Chloridoid/PACCAD
Uniola paniculata / Chloridoid/PACCAD
Zoysia / Chloridoid/PACCAD
53
Spartina / Chloridoid/PACCAD
42
Sporobolus / Chloridoid/PACCADDistichlis spicata / Chloridoid/PACCAD

Pappophorum g. / Chloridoid/PACCAD
Micraira / i.s.
Eriachne / i.s.
Gynerium sagittatum / i.s.
Oryza sativa / Ehrhartoid/BEP
67
Leersia / Ehrhartoid/BEP

51

63
49

Danthonia / Danthonioid/PACCAD

Merxmuellera macowanii / Danthonioid/PACCAD
Molinia g. / Arundinoid/PACCAD

56
43

Karoochloa purpurea / Danthonioid/PACCAD
Austrodanthonia laevis / Danthonioid/PACCAD

Ehrharta / Ehrhartoid/BEP
Streptogyna americana / i.s.
98
50

68

Lygeum spartum / Pooid/BEP
Nardus stricta / Pooid/BEP
Melica / Pooid/BEP

Glyceria / Pooid/BEP
Bromus / Pooid/BEP
Triticum g. / Pooid/BEP
52
Avena g. / Pooid/BEP
53
Brachypodium / Pooid/BEP
Ampelodesmos mauritanica / Pooid/BEP
52
Piptatherum / Pooid/BEP
Phaenosperma globosum / Pooid/BEP
Brachyelytrum erectum / Pooid/BEP
Anisopogon avenaceus / Pooid/BEP
Stipa / Pooid/BEP
Nassella / Pooid/BEP
Diarrhena / Pooid/BEP
77

88
50

49
49

Pseudosasa g. / Bambusoid/BEP
Chusquea / Bambusoid/BEP
Eremitis sp. / Bambusoid/BEP
93
Pariana radiciflora / Bambusoid/BEP
50
Lithachne / Bambusoid/BEP
61 Olyra / Bambusoid/BEP
46
Buergersiochloa bambusoides / Bambusoid/BEP

Joinvillea
ElegiaBaloskion tetraphyllum
0.05

Figure 7.8.: Total Evidence Poaceae phylogeny of the 66 taxa from the GPWG dataset.
For the analysis HKY model has been assumed, missing parameters havebeen estimated from the joint datasets. The numbers at the internal nodes
represent the percentage of the corresponding congruent splits.

80

7. Phylogenetic Trees from Multiple Genesets with Missing Data

Danthonia / Danthonioid/PACCAD
Karoochloa purpurea / Danthonioid/PACCAD
Austrodanthonia laevis / Danthonioid/PACCAD
Merxmuellera macowanii / Danthonioid/PACCAD
Aristida / Aristidoid/PACCAD
Stipagrostis zeyheri / Aristidoid/PACCAD
Miscanthus g. / Panicoid/PACCAD
Zea mays / Panicoid/PACCAD
Pennisetum / Panicoid/PACCAD
Panicum / Panicoid/PACCAD
Gynerium sagittatum / i.s.
Danthoniopsis / Panicoid/PACCAD
Zeugites pittieri / Centothecoid/PACCAD
Thysanolaena maxima / Centothecoid/PACCAD
Chasmanthium / Centothecoid/PACCAD
Molinia g. / Arundinoid/PACCAD
Phragmites australis / Arundinoid/PACCAD
Amphipogon / Arundinoid/PACCAD
Arundo donax / Arundinoid/PACCAD
Uniola paniculata / Chloridoid/PACCAD
Eragrostis / Chloridoid/PACCAD
Spartina / Chloridoid/PACCAD
Zoysia / Chloridoid/PACCAD
Sporobolus / Chloridoid/PACCADDistichlis spicata / Chloridoid/PACCAD
Pappophorum g. / Chloridoid/PACCAD
Merxmuellera rangei / Chloridoid/PACCAD
Centropodia glauca / Chloridoid/PACCAD
Micraira / i.s.
Eriachne / i.s.
Triticum g. / Pooid/BEP
Bromus / Pooid/BEP
Avena g. / Pooid/BEP
Brachypodium / Pooid/BEP
Diarrhena / Pooid/BEP
Stipa / Pooid/BEP
Nassella / Pooid/BEP
Ampelodesmos mauritanica / Pooid/BEP
Piptatherum / Pooid/BEP
Melica / Pooid/BEP
Glyceria / Pooid/BEP
Anisopogon avenaceus / Pooid/BEP
Phaenosperma globosum / Pooid/BEP
Nardus stricta / Pooid/BEP
Lygeum spartum / Pooid/BEP
Brachyelytrum erectum / Pooid/BEP
Pariana radiciflora / Bambusoid/BEP
Eremitis sp. / Bambusoid/BEP
Buergersiochloa bambusoides / Bambusoid/BEP
Lithachne / Bambusoid/BEP
Olyra / Bambusoid/BEP
Chusquea / Bambusoid/BEP
Pseudosasa g. / Bambusoid/BEP
Ehrharta / Ehrhartoid/BEP
Leersia / Ehrhartoid/BEP
Oryza sativa / Ehrhartoid/BEP
Streptogyna americana / i.s.
Guaduella marantifolia / Puelioid
Puelia ciliata / Puelioid
Pharus / Pharoid
Streptochaeta / Anomochlooid
Anomochloa marantoideai / Anomochlooid
Joinvillea
Baloskion tetraphyllum
ElegiaFlagellaria indica

Figure 7.9.: MRP phylogeny of the 66 taxa from the Poaceae dataset based on ML trees
and the Baum/Ragan coding scheme (MRP-BR).

7.4. The Phylogeny of the Grasses

81

Pappophorum g. / Chloridoid/PACCADDistichlis spicata / Chloridoid/PACCAD
Spartina / Chloridoid/PACCAD
Zoysia / Chloridoid/PACCAD
Sporobolus / Chloridoid/PACCAD
Uniola paniculata / Chloridoid/PACCAD
Eragrostis / Chloridoid/PACCAD
Merxmuellera rangei / Chloridoid/PACCAD
Eriachne / i.s.
Micraira / i.s.
Centropodia glauca / Chloridoid/PACCAD
Aristida / Aristidoid/PACCAD
Stipagrostis zeyheri / Aristidoid/PACCAD
Merxmuellera macowanii / Danthonioid/PACCAD
Danthonia / Danthonioid/PACCAD
Karoochloa purpurea / Danthonioid/PACCAD
Austrodanthonia laevis / Danthonioid/PACCAD
Arundo donax / Arundinoid/PACCAD
Amphipogon / Arundinoid/PACCAD
Gynerium sagittatum / i.s.
Chasmanthium / Centothecoid/PACCAD
Zeugites pittieri / Centothecoid/PACCAD
Thysanolaena maxima / Centothecoid/PACCAD
Danthoniopsis / Panicoid/PACCAD
Panicum / Panicoid/PACCAD
Pennisetum / Panicoid/PACCAD
Miscanthus g. / Panicoid/PACCAD
Zea mays / Panicoid/PACCAD
Phragmites australis / Arundinoid/PACCAD
Molinia g. / Arundinoid/PACCAD
Leersia / Ehrhartoid/BEP
Streptogyna americana / i.s.
Oryza sativa / Ehrhartoid/BEP
Ehrharta / Ehrhartoid/BEP
Chusquea / Bambusoid/BEP
Pseudosasa g. / Bambusoid/BEP
Buergersiochloa bambusoides / Bambusoid/BEP
Eremitis sp. / Bambusoid/BEP
Pariana radiciflora / Bambusoid/BEP
Lithachne / Bambusoid/BEP
Olyra / Bambusoid/BEP
Nardus stricta / Pooid/BEP
Lygeum spartum / Pooid/BEP
Triticum g. / Pooid/BEP
Bromus / Pooid/BEP
Avena g. / Pooid/BEP
Brachypodium / Pooid/BEP
Diarrhena / Pooid/BEP
Stipa / Pooid/BEP
Anisopogon avenaceus / Pooid/BEP
Piptatherum / Pooid/BEP
Ampelodesmos mauritanica / Pooid/BEP
Nassella / Pooid/BEP
Glyceria / Pooid/BEP
Melica / Pooid/BEP
Phaenosperma globosum / Pooid/BEP
Brachyelytrum erectum / Pooid/BEP
Puelia ciliata / Puelioid
Guaduella marantifolia / Puelioid
Pharus / Pharoid
Streptochaeta / Anomochlooid
Anomochloa marantoideai / Anomochlooid
Joinvillea
ElegiaBaloskion tetraphyllum
Flagellaria indica

Figure 7.10.: MRP phylogeny of the 66 taxa from the Poaceae dataset based on ML
trees and Purvis’ coding scheme (MRP-Pu).

82

7. Phylogenetic Trees from Multiple Genesets with Missing Data

PACCAD clade is again monophyletic. The PACCAD subfamilies Chloridoids, Aristidoids, and Danthonioids could be established where the latter two form a monophyletic
group. The branching pattern of Arundinoids, Panicoids, and Centothecoids remains
unresolved. The BEP clade consisting of the monophyletic Bambusoids, Erhartoids,
Pooids could be recovered whereas their branching order remained unresolved.
Using Purvis’ encoding one most parsimonious tree was reconstructed (MRP-Pu,
Fig. 7.10), which showed a by and large linearized topology. The PACCAD clade isrecovered in which the Chloridoids, Centothecoids, Aristidoids, and Danthonioids remain monophyletic. The members of the Arundinoids and Panicoids are placed on thetwo main PACCAD branches in a linearized order not providing separate subtrees. The
taxa of the BEP clade subfamilies are inserted along the main branch of grass subtree.
Only the Bambusoids subtree could be recovered.

7.4.5.3. MinCut Supertrees
The MinCut Supertree produced a long caterpillar tree (Fig. 7.11) that only maintains the Danthonioid subfamily at the top and a monophyletic group for the two earlybranching Anomochloids. In the linearized tree, the early branching taxa are separated
from the outgroup taxa and from the grass subtree by edges above and below.
The modifications suggested by Page (2002) in the ModMinCut Supertree algorithm lead to a biologically reasonable tree (Fig. 7.12). Not only the PACCAD clade
and the BEP clade are recovered, but also nearly all subfamilies could be established in
the highly resolved tree. In the PACCAD clade, the Panicoids are joined by the Centothecoids which do not show a monophyletic group. Then the other subfamily subtrees
join in the order Arundinoids, Danthonioids, Chloridoids, and Aristidoids. Within the
BEP clade Bambusoids and Pooids are sister groups followed by the Erhartoids. The
early branching taxa branch off from the root in the order Anomochloids, Pharus, andPuelioids.

7.4.5.4. Combined Quartet Tree
The reconstructed phylogeny of the grasses (Poaceae) based on the combined quartets is given in Fig. 7.13. Rooted by the outgroup taxa (Flagellaria, Baloskion, Elegia, and Joinvillea), the so-called early branching taxa Puelioids, Anomochloids, andPharus are placed in a multifurcation at the root of the grass subtree. The Puelioids
and Anomochloids form monophyletic groups. Within the grass subtree a group of six
subfamilies establishes the so-called PACCAD clade. Three of these subfamilies of the
PACCAD clade, the Aristidoids, Arundinoids, and Danthonioids, could be reconstructed,
while from the Panicoids, Centothecoids, Chloridoids only parts could be established asa monophyletic group.

7.4. The Phylogeny of the Grasses

83

c0
c2
c5
c7

c0
c0
c0
c0

Karoochloa purpurea / Danthonioid/PACCAD
Austrodanthonia laevis / Danthonioid/PACCAD
Danthonia / Danthonioid/PACCAD

Merxmuellera macowanii / Danthonioid/PACCAD

Merxmuellera rangei / Chloridoid/PACCAD

Eragrostis / Chloridoid/PACCAD

c9
Centropodia glauca / Chloridoid/PACCAD
c13
Arundo donax / Arundinoid/PACCAD
c16
Phragmites australis / Arundinoid/PACCAD
Molinia g. / Arundinoid/PACCAD
c16
Amphipogon / Arundinoid/PACCAD
c25
Miscanthus g. / Panicoid/PACCAD
Zea mays / Panicoid/PACCAD
c0
Pennisetum / Panicoid/PACCAD
c22
Gynerium sagittatum / i.s.
c24
Panicum / Panicoid/PACCAD
c31
Thysanolaena maxima / Centothecoid/PACCAD
c35
Chasmanthium / Centothecoid/PACCAD
Spartina / Chloridoid/PACCAD
c32
Micraira / i.s.
c38
Aristida / Aristidoid/PACCAD
c41
Stipagrostis zeyheri / Aristidoid/PACCAD
c24
Triticum g. / Pooid/BEP
c28
Bromus / Pooid/BEP
Avena g. / Pooid/BEP
c30
Stipa / Pooid/BEP
c44
Oryza sativa / Ehrhartoid/BEP
Lygeum spartum / Pooid/BEP
c48
Glyceria / Pooid/BEP
Melica / Pooid/BEP
c55
Ehrharta / Ehrhartoid/BEP
Nardus stricta / Pooid/BEP
Anisopogon avenaceus / Pooid/BEP
c56
Chusquea / Bambusoid/BEP
Lithachne / Bambusoid/BEP
c64
Sporobolus / Chloridoid/PACCAD
Nassella / Pooid/BEP
Diarrhena / Pooid/BEP
Brachypodium / Pooid/BEP
c62
Pappophorum g. / Chloridoid/PACCAD
c55
Leersia / Ehrhartoid/BEP
c58
Brachyelytrum erectum / Pooid/BEP
Piptatherum / Pooid/BEP
Ampelodesmos mauritanica / Pooid/BEP
c51
Pseudosasa g. / Bambusoid/BEP
Pariana radiciflora / Bambusoid/BEP
Eremitis sp. / Bambusoid/BEP
c50
Olyra / Bambusoid/BEP
c45
Danthoniopsis / Panicoid/PACCAD
Streptogyna americana / i.s.
Buergersiochloa bambusoides / Bambusoid/BEP
c26
Distichlis spicata / Chloridoid/PACCAD
Zoysia / Chloridoid/PACCAD
Uniola paniculata / Chloridoid/PACCAD
Zeugites pittieri / Centothecoid/PACCAD
Phaenosperma globosum / Pooid/BEP
c24
Puelia ciliata / Puelioid
Guaduella marantifolia / Puelioid
Pharus / Pharoid
Anomochloa marantoideai / Anomochlooid
Streptochaeta / Anomochlooid

Joinvillea

ElegiaBaloskion tetraphyllum
Flagellaria indica

Figure 7.11.: MinCut supertree of the 66 taxa from the Poaceae dataset based on the
8 separate ML trees. The numbers at the inner nodes present the numberof edges in the minimal cut set.

84

7. Phylogenetic Trees from Multiple Genesets with Missing Data

c0
c0
c0
hc1
c0
hc5
hc1

Zoysia / Chloridoid/PACCAD
Sporobolus / Chloridoid/PACCAD
Spartina / Chloridoid/PACCADDistichlis spicata / Chloridoid/PACCAD

Pappophorum g. / Chloridoid/PACCAD
Eragrostis / Chloridoid/PACCAD
Uniola paniculata / Chloridoid/PACCAD

Micraira / i.s.

Centropodia glauca / Chloridoid/PACCAD
Merxmuellera rangei / Chloridoid/PACCAD
c0
hc1

hc1

Karoochloa purpurea / Danthonioid/PACCAD
Austrodanthonia laevis / Danthonioid/PACCAD
Danthonia / Danthonioid/PACCAD

Merxmuellera macowanii / Danthonioid/PACCAD

hc4
hc3
hc1

Molinia g. / Arundinoid/PACCAD
Phragmites australis / Arundinoid/PACCAD
Arundo donax / Arundinoid/PACCAD
Amphipogon / Arundinoid/PACCAD
Chasmanthium / Centothecoid/PACCAD
c0
hc1

Thysanolaena maxima / Centothecoid/PACCAD
Danthoniopsis / Panicoid/PACCAD
Zeugites pittieri / Centothecoid/PACCAD
Miscanthus g. / Panicoid/PACCAD
Zea mays / Panicoid/PACCAD

c0
c0
hc2

Panicum / Panicoid/PACCAD
Pennisetum / Panicoid/PACCAD

Gynerium sagittatum / i.s.

Aristida / Aristidoid/PACCAD
Stipagrostis zeyheri / Aristidoid/PACCAD

hc1

c0
c0
hc2

Oryza sativa / Ehrhartoid/BEP
Leersia / Ehrhartoid/BEP
Ehrharta / Ehrhartoid/BEP

Streptogyna americana / i.s.
Pseudosasa g. / Bambusoid/BEP
Chusquea / Bambusoid/BEP

hc1
hc2

hc1
hc1

Buergersiochloa bambusoides / Bambusoid/BEP
Lithachne / Bambusoid/BEP
Olyra / Bambusoid/BEP
Eremitis sp. / Bambusoid/BEP
Pariana radiciflora / Bambusoid/BEP

Brachyelytrum erectum / Pooid/BEP
Lygeum spartum / Pooid/BEP
Nardus stricta / Pooid/BEP

hc1c0

Phaenosperma globosum / Pooid/BEP
hc1

Anisopogon avenaceus / Pooid/BEP
Ampelodesmos mauritanica / Pooid/BEP
Melica / Pooid/BEP
Glyceria / Pooid/BEP

c0
hc1
hc1
hc2

Stipa / Pooid/BEP
Nassella / Pooid/BEP

Piptatherum / Pooid/BEP
Diarrhena / Pooid/BEP
c0

Brachypodium / Pooid/BEP
c0

hc1
c0
c0
c0
c0

Bromus / Pooid/BEP
Triticum g. / Pooid/BEP
Avena g. / Pooid/BEP

Guaduella marantifolia / Puelioid
Puelia ciliata / Puelioid

Pharus / Pharoid
Anomochloa marantoideai / Anomochlooid
Streptochaeta / Anomochlooid

Joinvillea

ElegiaBaloskion tetraphyllum
Flagellaria indica

Figure 7.12.: ModMinCut supertree of the 66 taxa from the Poaceae dataset basedon the 8 separate ML trees. The numbers at the inner nodes present the
number of cuts or hierarchical cuts weighted as described by Page (2002).

7.4. The Phylogeny of the Grasses
Flagellaria indica
99
61
85

76

74

64

85

Guaduella marantifolia / Puelioid
Puelia ciliata / Puelioid
Anomochloa marantoideai / Anomochlooid
Streptochaeta / Anomochlooid
Aristida / Aristidoid/PACCAD
97
Stipagrostis zeyheri / Aristidoid/PACCAD
88
Molinia g. / Arundinoid/PACCAD
62
Phragmites australis / Arundinoid/PACCAD
51
Arundo donax / Arundinoid/PACCAD
47
Amphipogon / Arundinoid/PACCAD
Karoochloa purpurea / Danthonioid/PACCAD
80
Austrodanthonia laevis / Danthonioid/PACCAD
71
Danthonia / Danthonioid/PACCAD
66
Merxmuellera macowanii / Danthonioid/PACCAD
Chasmanthium / Centothecoid/PACCAD
52
Zeugites pittieri / Centothecoid/PACCAD
Thysanolaena maxima / Centothecoid/PACCAD
Pennisetum / Panicoid/PACCAD
58
Miscanthus g. / Panicoid/PACCAD
Danthoniopsis / Panicoid/PACCAD
83
Panicum / Panicoid/PACCAD
Zea mays / Panicoid/PACCADDistichlis spicata / Chloridoid/PACCAD
93
Pappophorum g. / Chloridoid/PACCAD
Spartina / Chloridoid/PACCAD
77
Sporobolus / Chloridoid/PACCAD
80
Zoysia / Chloridoid/PACCAD
Merxmuellera rangei / Chloridoid/PACCAD
78
Centropodia glauca / Chloridoid/PACCAD
Eragrostis / Chloridoid/PACCAD
60
Uniola paniculata / Chloridoid/PACCAD
53
Eriachne / i.s.
Micraira / i.s.
73
Gynerium sagittatum / i.s.
Oryza sativa / Ehrhartoid/BEP
96
Leersia / Ehrhartoid/BEP
88
Ehrharta / Ehrhartoid/BEP
Lygeum spartum / Pooid/BEP
99
Nardus stricta / Pooid/BEP
73
Bromus / Pooid/BEP
71
Triticum g. / Pooid/BEP
79
Brachypodium / Pooid/BEP
Diarrhena / Pooid/BEP
Avena g. / Pooid/BEP
65
Phaenosperma globosum / Pooid/BEP
96
Anisopogon avenaceus / Pooid/BEP
45
Ampelodesmos mauritanica / Pooid/BEP
55
Stipa / Pooid/BEP
60
Nassella / Pooid/BEP
50
Piptatherum / Pooid/BEP
Brachyelytrum erectum / Pooid/BEP
61
Glyceria / Pooid/BEP
55
Melica / Pooid/BEP
Eremitis sp. / Bambusoid/BEP
98
Pariana radiciflora / Bambusoid/BEP
91
Lithachne / Bambusoid/BEP
48
Olyra / Bambusoid/BEP
73
Buergersiochloa bambusoides / Bambusoid/BEP
54
Pseudosasa g. / Bambusoid/BEP
64
Streptogyna americana / i.s.
67
Chusquea / Bambusoid/BEP

Pharus / Pharoid
Joinvillea
ElegiaBaloskion tetraphyllum

Figure 7.13.: Combined Quartet phylogeny of the 66 taxa from the Poaceae dataset. For
the analysis the HKY model has been assumed, missing parameters havebeen estimated separately from the according datasets. The numbers at
the internal nodes represent the percentage of the corresponding congruent
splits.

86

7. Phylogenetic Trees from Multiple Genesets with Missing Data

From the BEP clade all three subfamilies, Bambusoids, Erhartoids, and Pooids, could
be established. The BEP clade, however, could not be reconstructed with the combined
quartet method.
Another consensus method, the so-called extended majority consensus, from the CONSENSE program (PHYLIP 3.6a3, Felsenstein, 1993) has been applied by others (Ben-dor
et al., 1998; Ranwez and Gascuel, 2001) to construct a binary tree from the set of intermediate trees. The extended majority consensus is a greedy consensus method that
adds all splits in descending order of percent occurrence to build a fully resolved tree.
All splits that contradict already accepted splits are discarded. We call this approach
Mext .
Applying Mext to the set of 100,000 intermediate trees produces a fully resolved tree
as anticipated (Fig. 7.14). In addition to the groups resolved by the Mrel consensus,
this tree reveals all PACCAD subfamilies, with the exception of the placement of Zea
mays. The Mext consensus tree presents a monophyletic BEP clade with Bambusoids
and Erhartoids being sister groups. The branching order of the early branching taxa is
resolved as Puelioids, Anomochloids, and finally Pharus.

7.4.6. Other Published Poaceae Phylogenies
Different versions of the Poaceae dataset have been analyzed by others (Grass Phylogeny
Working Group, 2001; Salamin et al., 2002). Since the Poaceae data is maintained by
the GPWG, it has grown over the last years. Therefore, the datasets used in publications differ in size (66 taxa in Grass Phylogeny Working Group, 2001 and 61 taxa in
Salamin et al., 2002) as possibly in the distribution of genes per taxon. Nevertheless, thepublished results will be discussed here, because the authors also used other methods toreconstruct combined phylogenies.
Three trees from the two publications (two total evidence trees and one supertree) aresketched in Fig. 7.15. The GPWG total evidence tree was reconstructed with maximum
parsimony from 8 datasets covering 66 taxa. In contrast to our study they used the
additional morphology and restriction data. Salamin et al. (2002) analyzed the 8 datasetswith a total of 61 taxa. Their total evidence tree was also constructed by MP, while the
MRP supertree (the best out of several different MRP trees) is based on a Baum/Ragan
encoding weighted by bootstrap values from the input trees.
All these trees show a similar overall topology having the early branching taxa at
the root of the grass subtree. All four support a monophyletic PACCAD clade. Also
Bambusoids, Erhartoids, and Pooids are shown as monophyletic subfamilies. Concerning
the order of the three subfamilies within the BEP clade, all three possible topologies
have been found. The GPWG total evidence tree groups together Bambusoids andEhrhartoids. The total evidence tree of Salamin et al. (2002) forms a Bambusoid/Pooid
group while their supertree supports a cluster of Ehrhartoids and Pooids. This shows thehigh amount of uncertainty in the data. The subfamily order within the PACCAD clade
remains an open question. All trees exhibit different branching orders, although most of7.4. The Phylogeny of the Grasses

87

66

Stipagrostis zeyheri / Aristidoid/PACCAD
Aristida / Aristidoid/PACCAD
Merxmuellera macowanii / Danthonioid/PACCAD
Danthonia / Danthonioid/PACCAD
71
Karoochloa purpurea / Danthonioid/PACCAD
80
Austrodanthonia laevis / Danthonioid/PACCAD

78

Merxmuellera rangei / Chloridoid/PACCAD
Centropodia glauca / Chloridoid/PACCAD

97

93

Distichlis spicata / Chloridoid/PACCAD
Pappophorum g. / Chloridoid/PACCAD

Uniola paniculata / Chloridoid/PACCAD
Eragrostis / Chloridoid/PACCAD
Eriachne / i.s.
Spartina / Chloridoid/PACCAD
77
Sporobolus / Chloridoid/PACCAD
Zoysia / Chloridoid/PACCAD
60

53

80
62
51

Molinia g. / Arundinoid/PACCAD
Phragmites australis / Arundinoid/PACCAD

Arundo donax / Arundinoid/PACCAD
Amphipogon / Arundinoid/PACCAD
Zea mays / Panicoid/PACCAD
Gynerium sagittatum / i.s.
73
Micraira / i.s.
Thysanolaena maxima / Centothecoid/PACCAD
Chasmanthium / Centothecoid/PACCAD
52
Zeugites pittieri / Centothecoid/PACCAD
Miscanthus g. / Panicoid/PACCAD
58
Pennisetum / Panicoid/PACCAD
Danthoniopsis / Panicoid/PACCAD
83
Panicum / Panicoid/PACCAD
Piptatherum / Pooid/BEP
50
Nassella / Pooid/BEP
60
Stipa / Pooid/BEP
55
Ampelodesmos mauritanica / Pooid/BEP
Anisopogon avenaceus / Pooid/BEP
96
Phaenosperma globosum / Pooid/BEP
45
Melica / Pooid/BEP
55
Glyceria / Pooid/BEP
61
Brachyelytrum erectum / Pooid/BEP
65
Brachypodium / Pooid/BEP
Bromus / Pooid/BEP
71
Triticum g. / Pooid/BEP
Diarrhena / Pooid/BEP
Avena g. / Pooid/BEP
Lygeum spartum / Pooid/BEP
99
Nardus stricta / Pooid/BEP
47

88

76

73

Lithachne / Bambusoid/BEP
Olyra / Bambusoid/BEP
73
Buergersiochloa bambusoides / Bambusoid/BEP
91
Eremitis sp. / Bambusoid/BEP
98
Pariana radiciflora / Bambusoid/BEP
Chusquea / Bambusoid/BEP
67
Streptogyna americana / i.s.
64
Pseudosasa g. / Bambusoid/BEP
Ehrharta / Ehrhartoid/BEP
Oryza sativa / Ehrhartoid/BEP
96
Leersia / Ehrhartoid/BEP
48

54

88

85
61
99
74

Pharus / Pharoid
Anomochloa marantoideai / Anomochlooid

Streptochaeta / Anomochlooid
Puelia ciliata / Puelioid
Guaduella marantifolia / Puelioid

Joinvillea
Elegia64
Baloskion tetraphyllum
Flagellaria indica

Figure 7.14.: Extended majority consensus of all intermediate trees from the combined
reconstruction. The numbers at the internal nodes represent the percentageof the corresponding congruent splits.

88

7. Phylogenetic Trees from Multiple Genesets with Missing Data
Outgroup
Outgroup
Outgroup

Early branching taxa
Early branching taxa

97
36

100
100

100

61

Pooideae

Erhartoideae

87

100

47

Arundinoideae

86

Chloridoideae
Centothecoideae

68

Gynerium (i.s.)

68
65

Panicoideae

Erhartoideae
Pooideae

Erhartoideae

Aristidoideae

Chloridoideae
100

Danthonioideae

77

Streptogyna (i.s.)

Streptogyna (i.s.)

100

Aristidoideae
100
98

32

87

a

98

Pooideae
Micraira (i.s.)

7125

Bambusoideae

Streptogyna (i.s.)

40

100

98
96−100Bambusoideae

100

Bambusoideae

100

Early branching taxa
100

52

Micraira (i.s.)
Arundinoideae

100

Aristidoideae

Panicoideae+Gynerium (i.s.)+Centothecoideae
Arundinoideae

Danthonioideae
76b

Chloridoideae

Centothecoideae+ Gynerium (i.s.)
Panicoideae

DanthoniacMicraira (i.s.)

Figure 7.15.: Poaceae phylogenies from other publications: (a) total evidence tree from
Grass Phylogeny Working Group (2001); (b) total evidence tree, and (c)supertree from Salamin et al. (2002). The numbers represent bootstrap
values.
the subfamilies are displayed as monophyletic. In the GPWG and Salamin et al.’s total
evidence trees, Aristidoids and Danthonioids form sister groups. The same holds forCentothecoids and Panicoids. The topology of the PACCAD subtree in Salamin et al.’s
supertree differs strongly from the others.

7.5. Discussion
7.5.1. Problems of Dataset-Combining Methods
There is a long ongoing debate on the question whether low level or high level combination should be preferred (e.g. de Queiroz et al., 1995, and references therein). A mainargument against total evidence methods is the problem of how to choose an appropriate
model of evolution for the concatenated dataset. It is well-known that different areas of
the genome develop under different evolutionary constraints. It is hard if not impossible
to choose one appropriate model for different sequences coding for proteins, structuralor functional RNAs, and ’function-less’ parts of the genome. Missing data also causes
problems because many methods discard alignment columns containing gaps. In thiscase only few or even no alignment columns remain for the analysis.
Since consensus and supertree methods reconstruct trees separately for each sequencealignment, they can easily use different evolutionary models matching the different constraints of the genesets. Supertree methods have been developed to handle the problem of
missing data. For the reconstruction of an overall tree, consensus and supertree methods
are using meta data like Matrix Representation derived from the input trees, or they use
the trees directly. The sequence data is usually not used in this process. Although this7.5. Discussion

89

is an advantage when input trees obtained from publications are combined, it discards
valuable information when the underlying sequence data is available.
It has been stated by Bininda-Emonds (2002, 275) that ”the inherent loss of information from using the source trees is a necessary trade-off to be able to combine all possible
sources of phylogenetic information.” However, the medium level method proposed inthis chapter combines the datasets at the quartet level guided by the likelihood of underlying data and, thus, tries to minimize the loss of information in the combination process.
Furthermore, the method offers the possibility to use appropriate models of evolution forthe different sources of data. Yang (1996) has shown that combined analysis of datasets
is possible even if the parameters as well as the branch lengths differ among the datasets.
In general, it should also be possible to combine likelihoods of different types of data like
DNA, amino acid, and restriction data. A prerequisite for combining such different data
sources is the availability of applicable ML models to compute the quartet topologies.
The necessity to adequately normalize the substitution models for different types of dataremains to be elucidated.

7.5.2. Comparison of the Tree Topologies
The maximum likelihood framework is the usual way to test whether tree topologies are
significantly different and which tree might be supported best by the data. A numberof possible tests have been suggested (for review see Goldman et al., 2000).
Yet such tests require a model of evolution as well as a dataset that enables a comparison of topologies. In our example only overlapping incomplete genesets are available.
Substantial variation has been found among the parameter sets that have been estimated from the genesets (cf. section 7.4.3). Consequently, as discussed above, it might
be impossible to choose one parameter set that is adequate for the overall dataset. Moreover, for testing evolutionary hypotheses it is highly questionable to use a concatenateddataset with such a large amount of missing data. The concatenated GPWG dataset
contains 37.85% gaps and missing characters, while the amount of ambiguous characters
per sequence range from 5.42% (Zea mays) up to 76.41% for Eriachne.
Since this framework seems not applicable (cf. also Novacek, 2001), the comparison of
the presented trees which are based on various optimality criteria relies on the taxonomic
classification proposed in a joint effort by a group of experts, the Grass Phylogeny
Working Group (2001).
The essence of the classification is the following. The PACCAD clade is monophyletic
consisting of the six monophyletic subfamilies Aristidoids, Danthonioids, Chloridoids,
Centothecoids, Arundinoids, Panicoids as well as the three species incertae sedis (i.s.1 ).
The species i.s. are Eriachne, Micraira, and Gynerium. Three grass subfamilies, namely,
Bambusoids, Erhartoids, and Pooids are reported to form a monophyletic group for different datasets, the so-called BEP clade. Finally, three early branching subfamilies, the1latin for ’of uncertain placement’90

7. Phylogenetic Trees from Multiple Genesets with Missing Data

Anomochloids, Pharoids, and Puelioids were described. This comprehensive subfamilyclassification consists of eleven previously published subfamilies and the new Danthonioid
subfamily (see Grass Phylogeny Working Group, 2001). This taxonomic classification isbased on the review of a large number of publications about Poaceae as well as a Total
evidence tree and the separate genetrees from the GPWG dataset.
Examining the reconstructed trees shows that most methods can more or less resolve
the taxonomic structure of the grasses. Taxa from the unresolved subfamilies are generally placed in multifurcations which do not contradict the taxonomic classification.
Only MinCut Supertree and MRP-Pu distribute members of unrevealed subfamilies
throughout the constructed tree. The methods which where able to reconstruct best thegrass subfamilies, PACCAD, and BEP clade were ModMinCut Supertree, the Mext
consensus from the combined quartet intermediate trees, and the two MP-based total
evidence trees (Fig. 7.15).
Although the Mext consensus gave a biologically reasonable result, its application
should not generally be recommended. It has been suggested by Kluge and Kolf (1993)
that bold, i.e., completely resolved, hypothesis are to be preferred over ’safe’ consensus
trees, since the safest tree would be the unresolved one. In order to achieve a completely
resolved tree, however, the Mext consensus accepts splits which would have been outvoted
by more abundant splits if the latter would not have been discarded due to earlier
contradictions. Using the Mrel consensus instead avoids this problem, since it does not
resolve the topology beyond the first incongruence. This agrees more with de Queiroz
et al. (1995, 664), stating that ”one might want a tree in which all clades have a certain
level of support.”

7.5.3. Conclusion
We have presented a new quartet-based medium level method to combine multiple sequence datasets with missing data for phylogeny reconstruction. The results obtained
from the GPWG Poaceae dataset show that this method performs comparable to othermethods for phylogeny reconstruction from overlapping datasets with missing data. Although the resulting tree was not fully resolved, it did not contradict the taxonomic
classification by the Grass Phylogeny Working Group (2001).
The modifications to the MinCut Supertree suggested by Page (2002, ModMinCut Supertree) have shown to substantially increase the result (cf. 7.4.5.3). Yet,none of the method to analyze incomplete phylogenetic data presented here has provento be the final choice. Facing the problem that it is unlikely to have complete molecularand morphological datasets in the near future, methods combining overlapping incomplete datasets remain a valuable tool to reconstruct evolution from multiple datasets(cf. also Bininda-Emonds, 2002). The development and improvement of methods forthe combined analysis of multiple dataset should be encouraged. Tree reconstruction
using combined quartets provides a new tool for such work.

7.5. Discussion

91

The judgment of the efficiency of dataset combining methods highly relies on biological knowledge. Since the ’true’ tree is unknown in general for these complex types ofdatasets, we suggest to perform simulations to study the reliability of such methods.
Furthermore, the influence of composite taxa should be examined in more detail. Recently, Malia Jr. et al. (2003) presented results on misleading effects of composite taxa
on total evidence trees.
The medium level approach using combined quartets also offers a possibility to analyze heterogeneous datasets as follows. Fast evolving datasets might not be suitable toreconstruct deep phylogenies due to saturation, but might contain valuable informationof the relationship of closely related species. Such datasets can be decomposed into separate group-specific ’genesets’ which are then connected by datasets from slower evolving
genes.
In addition to tree reconstruction, the set of combined quartet tree likelihoods can also
be used for likelihood mapping analysis to visualize, for instance, the combined quartet
support for specified cluster relationships (cf. 3.2.1).

92

7. Phylogenetic Trees from Multiple Genesets with Missing Data

8. Summary
Efficient sequencing techniques lead to exponentially growing amounts of molecular sequence data resulting mainly in two kinds of large datasets: vertically (large number oftaxa) as well as horizontally large datasets (long or many sequences per taxon). The
eminent need for efficient methods to analyze such datasets motivated the development
and improvement of several approaches based on the Quartet Puzzling algorithm whichwere proposed in this thesis.
Accelerated Puzzling Step (chapter 4): Two efficient O(n4 ) algorithms, one splitbased and a recursive algorithm, were developed to accelerate the puzzling step since
it demands a major runtime portion of the tree reconstruction by Quartet Puzzling.
Runtime tests show that the new methods outperform the original as well as Berry’s
algorithm. The best performing algorithm is based on the recursive computation of theedge penalties. It requires about half the runtime of the original O(n5 ) algorithm.
Parallelized Quartet Puzzling (chapter 5): To reduce the overall waiting time forthe reconstruction of large trees, a parallel implementation was developed. For theparallel workflow we investigated several scheduling algorithms. At last, a GSS-basedscheduling algorithm, modified to our needs, was devised to produce an even workloadon heterogeneous workstation clusters which are the most abundant parallel platforms
in molecular biology. Benchmark tests of the parallel implementation showed convincing
scaling behavior for large datasets.
Phylogenetic Trees from Less Quartets (chapter 6): A second component of the
QP algorithm which demands a large part of computation time is the ML step. An
approach is suggested to reduce the number of required quartets. To that end thedataset is decomposed into overlapping subsets from which the sets of possible quartetsare assembled into an overall tree applying a modified puzzling step algorithm called
ModPUZZLE. Results show that the resolution one can gain depends on the overlap
between the subsets which can be determined by the number of quartets one is willing tocompute. The approach can be applied to reconstruct the rough tree topology of large
datasets.
Combined Quartet Method for Overlapping Datasets (chapter 7): We presented
a procedure to reconstruct overlapping datasets from different genes (genesets) with

93

94

8. Summary

missing data. In contrary to total evidence approaches our method combines the datasetson a so-called medium level allowing to use evolutionary parameters specific for eachgeneset. The combination of the data is performed on the quartet level guided by the
likelihoods computed from the sequences. Proceeding this way tries to minimize theloss of primary information before combining the genesets, a common criticism of thesupertree methods. Different methods to combine genesets were applied to the biologicaldataset of the grasses. With regard to the GPWG taxonomic classification of the grasses,
the combined quartet method performs comparable to the other methods.

A. Variables and Functions
e– taxa, on the left of edge ee– an edgee– taxa, on the right of edge ee(•,a)

r e– taxa, on the same side as theroot of edge e– edge ending at node a in arooted treee(a,b)– edge connecting nodes a and b

E– set of edgesEmin– edge set with minimum
penaltyr e– taxa, not on the same side as
the root of edge e∅– empty set

α– shape parameter of a
Γ-distribution

Emin (T ) – edge set with minimum
penalty in tree TA– character alphabetEa– a leaf labelA– alignment

E(. . .) – penalty induced by external
node . . .

Ac– column c of alignment A

EMP– evolutionary Markov process

Acj– character state in alignment
column c of sequence jF– front set of nodesg– a gene; a geneset index

ABS– atomic batch size (scheduling)G = (V, E) – graphb– a leaf label

Govl– overlap graphB(n)– number of possible binary
trees with n leaves

GST– glutathione S-transferasec– a leaf labelh– height of a nodehmax

– maximal height in a treeH– hypothesisi– general purpose variable;
insertion counterj– general purpose variable

C e (. . .) – cluster of edge . . .
C n (. . .) – cluster of node . . .
CLP (. . .) – sum of all penalty paths
beginning at any leaf incluster . . .

– external node penalty vectord– a leaf labelk– number of subsets/genesetsD– square distance`T– log-likelihood of TD

– distance matrix

L(T )– leaf set of tree T

Deg(. . .) – degree of node . . .

l– a node or leaf

εL(. . .) – likelihood of . . .

– an edge of the penalty path

95

96

Li (j)LTc
LA
T– likelihood of state j at node i– likelihood of T
– likelihood of T in alignment
column Acµ– calibration factor for
substitution matrix Qm– number of sites, sequencelengthM– evolutionary modelMl– l-majority consensus
Mext– extended majority consensus(PHYLIP)Mrel– relative majority consensus
MRCA(. . .) – most recent common
ancestor of . . .
ν– size of subsetsn– number of sequencesN– penalty neighborhood matrix
Na,b– entry of matrix NNT– number of tasks (scheduling)NP– number of processes(scheduling)NR– number of un-scheduled tasks(scheduling)o– overlap between subsetsΩ(. . .) – complexity measure, lowerboundO(. . .) – complexity measure, upperboundΦ(. . .) – penalty of edge . . .
π– frequency vector, stationary
distribution
πi– frequency of state ipi– Bayesian weight of topology iP– substitution probability
matrixPρσ– entry of P for substitution
ρ→σ

A. Variables and FunctionsPρσ (t) – probability of substitution
ρ → σ in time t
P r(. . .) – probability of . . .
q– quartetQ– set of all possible quartetsfrom mathcalSQi– set of all possible quartetsfrom mathcalSi

Q(X,Ti ) – set of all quartets containing
X and three sequences from Ti
Q– instantaneous rate matrix

Qij– entry of Q for substitutioni→jρ

– a character stateR– substitution rate matrixr– a node or leaf

R– root nodeRij

– entry of R for substitutioni→jσ– a character stateS– set of sequencesSi– subset of sequences

S1∩2– intersection of S1 and S2

si– ith sequenceST– graph from the input trees T(MinCut)

ST \ETM ax – graph from the input trees
T without ETM ax (MinCut)
SMS (. . .) – sum of all penalty paths
beginning at any leaf incluster . . . , while remaining
within the cluster
succ(. . .) – direct descendants of node...
Θ(. . .) – complexity measure, upperand lower bound97t– timeT (S)– tree containing leaf set St(a,b)– length of edge e(a,b) inevolutionary timeTn

– n-tree, tree with n leaves

Tnj– neighbor-joining treetMRCA – runtime (MRCA-based
puzzling step)u– a node or leafv

– a node or leafV– penalty propagation vector

torig– runtime (original puzzling
step)

trecur– runtime (recursive puzzling
step)

V

– set of nodes/verticesw– a node or leaf

tsplit– runtime (split-based puzzling
step)

wi

– generic weight of topology iT– a tree, tree topologyx– general purpose variable; index

T– set of trees

X– taxon to insert next

T |Si

– the trees from T with allleaves r 6∈ Si pruned(MinCut)Ξ– current status of a schedulery– general purpose variable; indexz– general purpose variable; index

B. Abbreviationsaa– amino acids

A– Adenine

ANSI

i.s.

– incertae sedis (lat.) – ofuncertain placement– American National Standard
Institute

ITS– internal transcribed spacer

JC69– Jukes-Cantor modelBR– Baum/Ragan encoding

JTT– Jones-Taylor-Thornton modelC– Cytosine

K2P– Kimura 2-Parameter model

COW– Cluster of Workstations

lslength – least square branch length
estimation

cpREV – reversible model for
chloroplast dataCPU

– Central Processing Unit,
Processor

LS len – least square branch length
estimation

DNA– Deoxyribonucleic Acid

MIMD – Multiple Instruction –Multiple DataED

– Evolutionary Distances

MinCut – Minimum Cut

EF– Elongation Factor

EMBL – European Molecular Biology
Laboratories

ML– Maximum Likelihood

mldistance – ML distance estimation
ML dist – ML distance estimation

EMP– evolutionary Markov process

ModMinCut – Modified Minimum Cut

G– GuanineMP– Maximum Parsimony

MPI– Message Passing Interface

GBSSI – Granule Bound Starch
Synthase I
GPWG – Grass Phylogeny Working
Group
GTR

– General Time Reversible
model

GSS

– Guided Self-Scheduling

HKY– Hasegawa-Kishino-Yano
model

MRCA – Most Recent Common
Ancestor
MRF– Matrix Representation using
Flipping

MRP– Matrix Representation using
Parsimony

MSA– Multiple Sequence Alignment

MST– Minimum Spanning Tree

HKY85 – Hasegawa-Kishino-Yano
model

mtREV – reversible model for
mitochondrial data

Indel

nt

98

– Insertions and Deletions– nucleotides

99

ndhF

– NADH dehydrogenase,subunit F

SGSS– Smooth Guided
Self-SchedulingPC– Personal Computer

SIMD

PE

– Processing Element– Single Instruction – Multiple
Data

phyB

– phytochrome B

SMP– Scalable Multi-Processor

Pu– Purvis’ encoding

QP– Quartet Puzzling

rbcL– ribulose 1,5-bisphosphatecarboxylase/oxygenase, large
subunitSPRNG – Scalable Parallel Random
Number Generator
SS– Self-Scheduling

ssu rRNA – small subunit rRNAT

– Thymine

RDP– Ribosomal Database Project

TN93– Tamura-Nei modelRNA– Ribonucleic acidts

– transitionrRNA– ribosomal ribonucleic acid

TSS

– Trapezoid Self-SchedulingrpoC– RNA polymerase II, β”
subunittv– transversion

U– Uracil– RNA polymerase II, β”
subunitUV– ultravioletVT– Variable Time model– Static Chunking (scheduling)

WAG– Whelan-and-Goldman model

rpoC2
SC

100

B. Abbreviations

Bibliography
Adachi, J. and Hasegawa, M. (1996a) Model of amino acid substitution in proteinsencoded by mitochondrial DNA. J. Mol. Evol., 42, 459–468.
Adachi, J. and Hasegawa, M. (1996b) MOLPHY Version 2.3 – Programs for MolecularPhylogenetics Based on Maximum Likelihood, vol. 28 of Computer Science Monographs. Institute of Statistical Mathematics, Minato-ku, Tokyo.
Adachi, J., Waddell, P. J., Martin, W. and Hasegawa, M. (2000) Plastid genome phylogeny and a model of amino acid substitution for proteins encoded by chloroplast
DNA. J. Mol. Evol., 50, 348–358.
Adams III, E. N. (1986) N-trees as nestings: Complexity, similarity, and consensus. J.
Classif., 3, 299–317.
Aho, A. V., Sagiv, Y., Szymanski, T. G. and Ullman, J. D. (1981) Inferring a tree
from lowest common ancestors with an application to the optimization of relational
expressions. SIAM J. Comput., 10, 405–421.
Ansorge, W., Zimmermann, J., Erfle, H., Hewitt, N., Rupp, T., Schwager, C., Sproat, B.,
Stegemann, J. and Voss, H. (1993) Sequencing reactions for ALF (EMBL) automated
DNA sequencer. Methods Mol. Biol., 23, 317–356.
Bandelt, H.-J. and Dress, A. (1986) Reconstructing the shape of a tree from observed
dissimilarity da ta. Adv. Appl. Math., 7, 309–343.
Barthelemy, J. P. and Guenoche, A. (1991) Trees and Proximity Representations. Interscience Series in Discrete Mathematics and Optimization, Wiley, New York.
Baum, B. R. (1992) Combining trees as a way of combining data sets for phylogenetic
inference, and the desirability of combining gene trees. Taxon, 41, 3–10.
Baxevanis, A. D. (2003) The molecular biology database collection: 2003 update. Nucleic
Acids. Res., 31, 1–12.
Baxevanis, A. D., Davison, D. B., Page, R. D. M., Stormo, G. and Stein, L. (eds.) (2002)
Current Protocols in Bioinformatics. Wiley and Sons, New York, USA.

101

102

Bibliography

Ben-dor, A., Chor, B., Graur, D., Ron, O. and Pelleg, D. (1998) Constructing phylogenies from quartets: Elucidation of eutherian superordinal relationships. J. Comput.
Biol., 5, 377–390.
Benson, D. A., Karsch-Mizrachi, I., Lipman, D. J., Ostell, J. and Wheeler, D. L. (2003)
GenBank. Nucleic Acids Res., 31, 23–27.
Bininda-Emonds, O. R. P. (2002) MRP supertree construction in the consensus setting.
In Bioconsensus: Proceedings of Tutorial and Workshop on Bioconsensus II, DIMACS
Series in Discrete Mathematics and Theoretical Computer Science, page N.N.
Blanchette, M., Bourque, G. and Sankoff, D. (1997) Breakpoint phylogenies. In Proceedings of the Genome Informatics Workshop VIII, pages 25–34, Universal Academy
Press.
Boeckmann, B., Bairoch, A., Apweiler, R., Blatter, M.-C., Estreicher, A., Gasteiger,
E., Martin, M. J., Michoud, K., O’Donovan, C., Phan, I., Pilbout, S. and Schneider,M. (2003) The SWISS-PROT protein knowledgebase and its supplement TrEMBL in
2003. Nucleic Acids Res., 31, 365–370.
Bremer, K. (1990) Combinable component consensus. Cladistics, 6, 369–372.
Bryant, D. and Steel, M. (1995) Extension operations on sets of leaf-labeled trees. Adv.
Appl. Math., 16, 425–453.
Chen, D., Diao, L., Eulenstein, O., Fernández-Baca, D. and Sanderson, M. J. (2003)
Flipping: A supertree construction method. DIMACS Series in Discrete Mathematics
and Theoretical Computer Sciences, 000, to appear.
Cormen, T. H., Leiserson, C. E., Rivest, R. L. and Stein, C. (2001) Introduction to
Algorithms. The MIT Press, Cambridge, Massachusetts, Second edn..
Cramer, J. S. (1989) Econometric Applications of Maximum Likelihood Methods. Cambridge University Press, Cambridge, UK.
Day, W. H. E. (1987) Computational complexity of inferring phylogenies from dissimilarity matrices. Bull. Math. Biol., 49, 461–467.
Day, W. H. E. and Sankoff, D. (1986) Computational complexity of inferring phylogenies
by compatibility. Syst. Zool., 35, 224–229.
Dayhoff, M. O., Schwartz, R. M. and Orcutt, B. C. (1978) A model of evolutionarychange in proteins. In Dayhoff, M. O. (ed.), Atlas of Protein Sequence Structure,
vol. 5, pages 345–352, National Biomedical Research Foundation, Washington DC.
Dempster, A. P., Laird, N. M. and Rubin, D. B. (1977) Maximum likelihood fromincomplete data via the em algorithm. J. Royal Statist. Soc., 39, 1–39.

Bibliography

103

Dress, A., von Haeseler, A. and Krüger, M. (1986) Reconstructing phylogenetic trees
using variants of the four-point condition. Studien zur Klassifikation, 17, 299–305.
Durbin, R., Eddy, S. R., Krogh, A. and Mitchison, G. (1998) Biological sequence analysis - Probabilistic models of proteins and nucleic acids. Cambridge University Press,
Cambridge.
Duret, L., Mouchiroud, D. and Gouy, M. (1994) HOVERGEN, a database of homologous
vertebrate genes. Nucleic Acids Res., 22, 2360–2365.
Edman, P. (1950) Method for determinatin of the amino acid sequence in peptides. Acta
Chem. Scand., 4, 283–290.
Edman, P. and Begg, G. (1967) A protein sequenator. Eur. J. Biochem., 1, 80–91.
Edmiston, E. and Wagner, R. A. (1987) Parallelization of the dynamic programmingalgorithm for comparison of sequences. In Proceedings of the 1987 International Conference on Parallel Processing, pages 78–80, Pennsylvania, Penn State Press.
El-Rewini, H., Lewis, T. G. and Ali, H. H. (1994) Task Scheduling in Parallel andDistributed Systems. Prentice-Hall, Englewood Cliffs, New Jersey.
Eliason, S. R. (ed.) (1993) Maximum Likelihood Estimation: Logic and Practice. Quantitative Applications in the Social Sciences, Sage Publications, Newbury Park.
Ewens, W. J. and Grant, G. R. (2001) Statistical Methods in Bioinformatics: An Introduction. Springer Verlag, New York, USA.
Felsenstein, J. (1978) The number of evolutionary trees. Syst. Zool., 27, 27–33.
Felsenstein, J. (1981) Evolutionary trees from DNA sequences: a maximum likelihood
approach. J. Mol. Evol., 17, 368–376.
Felsenstein, J. (1993) PHYLIP (Phylogeny Inference Package) version 3.5c. Departmentof Genetics, University of Washington, Seattle, Distributed by the author.
Fitch, W. M. (1981) A non-sequential method for constructing trees and hierarchicalclassifications. J. Mol. Evol., 18, 30–37.
Foulds, L. R. and Graham, R. L. (1982) The Steiner problem in phylogeny is NPcomplete. Adv. Appl. Math., 3, 43–49.
Goldman, N. (1990) Maximum likelihood inference of phylogenetic trees, with specialreference to a poisson process model of DNA substitution and to parsimony analyses.
Syst. Zool., 39, 345–361.
Goldman, N., Anderson, J. P. and Rodrigo, A. G. (2000) Likelihood-based tests of
topologies in phylogenetics. Syst. Biol., 49, 652–670.

104

Bibliography

Graham, R. L. and Foulds, L. R. (1982) Unlikelihood that minimal phylogenies for a
realistic biological study can be constructed in reasonable computational time. Math.
Biosci., 60, 133–142.
Grass Phylogeny Working Group (2001) Phylogeny and subfamilial classification of thegrasses (poaceae). Ann. Mo. Bot. Gard., 88, 373–457.
Graur, D. and Li, W.-H. (2000) Fundamentals of Molecular Evolution. Sinauer Associates, Sunderland, Massachusetts, Second edn..
Gropp, W., Huss-Lederman, S., Lumsdaine, A., Lusk, E., Nitzberg, B., Saphir, W. andSnir, M. (1998) MPI: The Complete Reference - The MPI Extensions, vol. 2. The MIT
Press, Cambridge, Massachusetts, Second edn..
Gu, X., Fu, Y.-X. and Li, W.-H. (1995) Maximum likelihood estimation of the heterogeneity of substitution rate among nucleotide sites. Mol. Biol. Evol., 12, 546–557.
Hagerup, T. (1997) Allocating independent tasks to parallel processors: An experimental
study. J. Parallel Distrib. Comput., 47, 185–197.
Hagstrom, R., Matsuda, H., Overbeek, R., Olsen, G. and Woese, C. (1992) Inferring relationships among microorganisms by using maximum likelihood. The Concurrent Supercomputing Consortium Annual Report FY1991-1992, California Institute of Technology, Pasadena, CA, USA.
Hasegawa, M., Kishino, H. and Yano, T.-A. (1985) Dating of the human–ape splitting
by a molecular clock of mitochondrial DNA. J. Mol. Evol., 22, 160–174.
Huang, X. (1989) A space-efficient parallel sequence comparison algorithm for a messagepassing multiprocessor. Int. J. Parallel Program., 18, 223–239.
Huelsenbeck, J. P. and Ronquist, F. (2001) MRBAYES: Bayesian inference of phylogenetic trees. Bioinformatics, 17, 754–755.
Huson, D. H., Nettles, S. M. and Warnow, T. J. (1999) Disk-covering, a fast-converging
method for phylogenetic reconstruction. J. Comput. Biol., 6, 369–386.
Jones, D. T., Taylor, W. R. and Thornton, J. M. (1992) The rapid generation of mutation
data matrices from protein sequences. Comput. Appl. Biosci., 8, 275–282.
Jukes, T. H. and Cantor, C. R. (1969) Evolution of protein molecules. In Munro, H. N.
(ed.), Mammalian Protein Metabolism, vol. 3, pages 21–123, Academic Press, NewYork.
Jülich, A. (1995) Implementations of BLAST for parallel computers. Comput. Appl.
Biosci., 11, 3–6.
Kimura, M. (1980) A simple method for estimating evolutionary rates of base substitutions through comparative studies of nucleotide sequences. J. Mol. Evol., 16, 111–120.

Bibliography

105

Kleinjung, J., Douglas, N. and Heringa, J. (2002) Parallelized multiple alignment. Bioinformatics, 18, 1270–1271.
Kluge, A. G. and Kolf, A. J. (1993) Cladistics: What’s in a word? Cladistics, 9, 183–199.
Knuth, D. E. (1997) Sorting and Searching, vol. 3 of The Art of Computer Programming.
Addison-Wesley, Reading, Massachusetts, Third edn..
Krause, A., Haas, S. A., Coward, E. and Vingron, M. (2002) SYSTERS, GeneNest,
SpliceNest: Exploring sequence space from genome to protein. Nucleic Acids Res.,30, 299–300.
Lander, E., Mesirov, J. P. and W., T. (1988) Protein sequence comparison on a data
parallel computer. In Proceedings of the 1988 International conference on Parallel
Processing, pages 257–263, Pennsylvania, Penn State Press.
Lewis, P. O. (1998) A genetic algorithm for maximum-likelihood phylogeny inference
using nucleotide sequence data. Mol. Biol. Evol., 15, 277–283.
Liébecq, C. (ed.) (1992) Biochemical Nomenclature and Related Documents. Portland
Press, London, Second edn..
Liu, J., Iba, H. and Ishizuka, M. (2001) Selecting informative genes with parallel genetic
algorithms in tissue classification. In Proceedings of the Genome Informatics Workshop
XII, vol. 12 of Genome Informatics, pages 14–23, Universal Academy Press.
Maidak, B. L., Cole, J. R., Lilburn, T. G., Parker Jr., C. T., Saxman, P. R., Farris,
R. J., Garrity, G. M., Olsen, G. J., Schmidt, T. M. and Tiedje, J. M. (2001) The
RDP-II (ribosomal database project). Nucleic Acids Res., 29, 173–174.
Malia Jr., M. J., Lipscomb, D. L. and Allard, M. W. (2003) The misleading effects ofcomposite taxa in supermatrices. Mol. Phylogenet. Evol., 27, 522–527.
Margush, T. and McMorris, F. R. (1981) Consensus n-trees. Bull. Math. Biol., 43, 239–
244.
McMorris, F. R. and Neumann, D. A. (1983) Consensus functions defined on trees. Math.
Soc. Sci, 4, 131–136.
Mendel, J. M. and Burrus, C. S. (1990) Maximum-Likelihood Deconvolution: A Journey
into Model-Based Signal Processing. Springer Verlag, New York.
Meyer, S. and von Haeseler, A. (2003) Identifying site-specific substitution. Mol. Biol.
Evol., 20, 182–189.
Meyer, S., Weiss, G. and von Haeseler, A. (1999) Pattern of nucleotide substitution andrate heterogeneity in the hypervariable regions I and II of human mtDNA. Genetics,
152, 1103–1110.

106

Bibliography

Morgenstern, B. (1999) DIALIGN 2: improvement of the segment-to-segment approach
to multiple sequence alignment. Bioinformatics, 15, 211–218.
Müller, T. and Vingron, M. (2000) Modeling amino acid replacement. J. Comput. Biol.,7, 761–776.
Nakaya, A., Yamamoto, K. and Yonezawa, A. (1995) RNA secondary structure prediction using highly parallel computers. Comput. Appl. Biosci., 11, 685–692.
Ness, S. R., Terpstra, W., Krzywinski, M., Marra, M. A. and Jones, S. J. M. (2002)
Assembly of fingerprint contigs: Parallelized FPC. Bioinformatics, 18, 484–485.
Novacek, M. J. (2001) Mammalian phylogeny: Genes and supertrees. Curr. Biol., 11,
R573–R575.
Olsen, G. J., Matsuda, H., Hagstrom, R. and Overbeek, R. (1994) fastDNAml: A toolfor construction of phylogenetic trees of DNA sequences using maximum likelihood.
Comput. Appl. Biosci., 10, 41–48.
Ottmann, T. and Widmayer, P. (2002) Algorithmen und Datenstrukturen. Spektrum
Akademischer Verlag, Heidelberg, Fourth edn..
Page, R. D. M. (2002) Modified mincut supertrees. In Proceedings of the 2nd Workshop
on Algorithms in Bioinformatics (WABI 2002), vol. 2452 of Lecture Notes in Computer
Science, pages 537–551, New York, Springer.
Page, R. D. M. and Holmes, E. C. (1998) Molecular Evolution: A Phylogenetic Approach.
Blackwell Science, Oxford.
Pedretti, K. T., Casavant, T. L., Braun, R. C., Scheetz, T. E., Birkett, C. L. and Roberts,
C. A. (1999) Three complementary approaches to parallelization of local BLAST service on workstation clusters. In Proceedings of the 5th International Conference on
Parallel Computing Technologies (PaCT-99), vol. 1662 of Lecture Notes in Computer
Science, pages 271–282, New York, Springer.
Van de Peer, Y., Baldauf, S. L., Doolittle, W. F. and Meyer, A. (2000a) An updated
and comprehensive rRNA phylogeny of (crown) eukaryotes based on rate-calibrated
evolutionary distances. J. Mol. Evol., 51, 565–576.
Van de Peer, Y., De Rijk, P., Wuyts, J., Winkelmans, T. and De Wachter, R. (2000b)
The European small subunit ribosomal RNA database. Nucleic Acids Res., 28, 175–
176.
Press, W. H., Flannery, B. P., Teukolsky, S. A. and Vetterling, W. T. (1988) Numerical
Recipes in C: The Art of Scientific Computing. Cambridge University Press, Cambridge.

Bibliography

107

Purvis, A. (1995) A composite estimate of primate phylogeny. Philos. Trans. R. Soc.
Lond. Ser. B, 348, 405–421.
de Queiroz, A., Donoghue, M. J. and Kim, J. (1995) Separate versus combined analysisof phylogenetic evidence. Annu. Rev. Ecol. Syst., 26, 657–681.
Ragan, M. A. (1992) Phylogenetic inference based on matrix representation of trees.
Mol. Phylogenet. Evol., 1, 53–58.
Rambaut, A. and Grassly, N. C. (1997) Seq-Gen: An application for the Monte
Carlo simulation of DNA sequence evolution along phylogenetic trees. Comput. Appl.
Biosci., 13, 235–238.
Ranwez, V. and Gascuel, O. (2001) Quartet-based phylogenetic inference: Improvements
and limits. Mol. Biol. Evol., 18, 1103–1116.
Robinson-Rechavi, M. and Graur, D. (2001) Usage optimization of unevenly sampled
data through the combination of quartet trees: An eutherian draft phylogeny basedon 640 nuclear and mitochondrial proteins. Isr. J. Zool., 47, 259–270.
Rodriguez, F., Oliver, J. L., Main, A. and Medina, J. R. (1990) The general stochastic
model of nucleotide substitution. J. theor. Biol., 142, 485–501.
Saitou, N. and Nei, M. (1987) The neighbor–joining method: A new method for reconstructing phylogenetic trees. Mol. Biol. Evol., 4, 406–425.
Salamin, N., Hodkinson, T. R. and Savolainen, V. (2002) Building supertrees: An empirical assessment using the grass family (poaceae). Syst. Biol., 51, 136–150.
Salter, L. A. and Pearl, D. K. (2001) Stochastic search strategy for estimation of maximum likelihood phylogenetic trees. Syst. Biol., 50, 7–17.
Sanderson, M. J., Purvis, A. and Henze, C. (1998) Phylogenetic supertrees: Assembling
the trees of life. TREE, 13, 105–109.
Sanger, F., Nicklen, S. and Coulsen, A. R. (1977) DNA sequencing with chain–terminating inhibitors. Proc. Natl. Acad. Sci. USA, 74, 5463–5467.
Sanger, F. and Thompson, E. O. P. (1953) The amino-acid sequence in the glycyl chainof insulin: 1. the identification of lower peptides from partial hydrolysates. Biochem.
J., 53, 353–365.
Sankoff, D. and Blanchette, M. (1998) Multiple genome rearrangement and breakpoint
phylogeny. J. Comput. Biol., 5, 555–570.
Sattath, S. and Tversky, A. (1977) Additive similarity trees. Psychometrika, 42, 319–
345.

108

Bibliography

Schmidt, H. A. (1996) Parallelisierung phylogenetischer Methoden zur Untersuchung der
Crown-Group-Radiation. Diplomarbeit, Universität zu Köln.
Schmidt, H. A. and von Haeseler, A. (2002) Quartet trees as a tool to reconstructlarge trees from sequences. In Jajuga, K., Sokolowski, A. and Bock, H.-H. (eds.), Data
Analysis, Classification, and Related Methods, Studies in Classification, Data Analysis,
and Knowledge Organization, pages 379–388, Springer, Heidelberg, New York.
Schmidt, H. A. and von Haeseler, A. (2003) Maximum likelihood analysis using TREEPUZZLE. In Baxevanis, A. D., Davison, D. B., Page, R. D. M., Stormo, G. and Stein,
L. (eds.), Current Protocols in Bioinformatics, pages 6.6.1–6.6.25, Wiley and Sons,
New York, USA.
Schmidt, H. A., Petzold, E., Vingron, M. and von Haeseler, A. (2003) Molecular phylogenetics: Parallelized parameter estimation and quartet puzzling. J. Parallel Distrib.
Comput., 63, in press.
Schmidt, H. A., Strimmer, K., Vingron, M. and von Haeseler, A. (2002) TREE-PUZZLE:
Maximum likelihood phylogenetic analysis using quartets and parallel computing.
Bioinformatics, 18, 502–504.
Semple, C. and Steel, M. (2000) A supertree method for rooted trees. Discr. Appl. Math.,
105, 147–158.
Shapiro, B. A., Jin Chu Wu, Bengali, D. and Potts, M. J. (2001) The massively parallel
genetic algorithm for RNA folding: MIMD implementation and population variation.
Bioinformatics, 17, 137–148.
Snir, M., Otto, S. W., Huss-Lederman, S., Walker, D. W. and Dongarra, J. (1998)
MPI: The Complete Reference - The MPI Core, vol. 1. The MIT Press, Cambridge,
Massachusetts, Second edn..
Stamatakis, A. P., Ludwig, T., Meier, H. and Wolf, M. J. (2002) AxML: A fast programfor sequential and parallel phylogenetic tree computations based on the maximum
likelihood method. In Proceedings of the 1st IEEE Computer Society Bioinformatics
Conference - CSB2002, pages 21–28, Palo Alto.
Stewart, C. A., Hart, D., Berry, D. K., Olsen, Gary J.and Wernert, E. A. and Fischer, W.
(2001) Parallel implementation and performance of fastDNAml - a program for maximum likelihood phylogenetic inference. In Proceedings of the International Conference
on High Performance Computing and Communications - SC2001, pages 191–201.
Stoesser, G., Baker, W., van den Broek, A., Garcia-Pastor, M., Kanz, C., Kulikova, T.,
Leinonen, R., Lin, Q., Lombard, V., Lopez, R., Mancuso, R., Nardone, F., Stoehr, P.,
Tuli, M. A., Tzouvara, K. and Vaughan, R. (2002) The EMBL nucleotide sequence
database: major new developments. Nucleic Acids Res., 31, 17–22.

Bibliography

109

Strimmer, K., Goldman, N. and von Haeseler, A. (1997) Bayesian probabilities andquartet puzzling. Mol. Biol. Evol., 14, 210–213.
Strimmer, K. and von Haeseler, A. (1996) Quartet puzzling: A quartet maximum–
likelihood method for reconstructing tree topologies. Mol. Biol. Evol., 13, 964–969.
Strimmer, K. and von Haeseler, A. (1997) Likelihood–mapping: A simple method tovisualize phylogenetic content of a sequence alignment. Proc. Natl. Acad. Sci. USA,94, 6815–6819.
Strimmer, K. S. (1997) Maximum Likelihood Methods in Molecular Phylogenetics. Ph.D.
thesis, Ludwig–Maximilians–Universität, München, Germany.
Swofford, D. L., Olsen, G. J., Waddell, P. J. and Hillis, D. M. (1996) Phylogeny reconstruction. In Hillis, D. M., Moritz, C. and Mable, B. K. (eds.), Molecular Systematics,
pages 407–514, Sinauer Associates, Sunderland, Massachusetts, Second edn..
Tamura, K. and Nei, M. (1993) Estimation of the number of nucleotide substitutionsin the control region of mitochondrial DNA in humans and chimpanzees. Mol. Biol.
Evol., 10, 512–526.
Tavare, S. (1986) Some probabilistic and statistical problems on the analysis of DNA
sequences. Lec. Math. Life Sci., 17, 57–86.
Thompson, J. D., Higgins, D. G. and Gibson, T. J. (1994) CLUSTAL W: Improving thesensitivity of progressive multiple sequence alignment through sequence weighting,positions–specific gap penalties and weight matrix choice. Nucleic Acids Res., 22,4673–4680.
Trelles, O. (2001) On the parallelisation of bioinformatics applications. Brief. Bioinform.,
2, 181–194.
Trelles, O., Ceron, C., Wang, H. C., Dopazo, J. and Carazo, J. M. (1998) Application
note: New phylogenetic venues opened by a novel implementation of the DNAml
algorithm. Bioinformatics, 14, 544–545.
Whelan, S. and Goldman, N. (2001) A general empirical model of protein evolution
derived from multiple protein families using a maximum likelihood approach. Mol.
Biol. Evol., 18, 691–699.
Yang, Z. (1996) Maximum-likelihood models for combined analyses of multiple sequence
data. J. Mol. Evol., 42, 587–596.
Yang, Z. and Kumar, S. (1994) Approximate methods for estimating the pattern of
nucleotide substitution and the variation of substitution rates among sites. Mol. Biol.
Evol., 13, 650–659.

110

Bibliography

Acknowledgments
First and foremost, I wish to thank my supervisor Arndt von Haeseler. I am gratefulfor his excellent advise and collaboration. Furthermore I want to thank him as well as
Ulrike and Jule Friedrichs for their extraordinary hospitality during my many visits in
Leipzig and Düsseldorf.
I want to thank Martin Vingron for introducing me to Arndt von Haeseler in thecontext of a joint project as well as for giving me the possibility the work in the pleasantand stimulating environment of his research groups at the DKFZ in Heidelberg and the
Max-Planck-Institute (MPI) for Molecular Genetics in Berlin.
Also I want to thank Bill Martin for accepting the task to read this thesis as a ’Korreferent’.
Special thanks go to Heike Hennig-Schmidt for taking the burden of reading loads oftext unrelated to her own field of research as well as to Antje Krause for patience torepeatedly listening to and discussing ideas.
I thank Korbinian Strimmer for a lot of help when getting used to his TREE-PUZZLE
source code and for communicating Vincent Berry’s MRCA-based algorithm.
I thank Ekkehard Petzold for his close collaboration on parameter estimation issues.
I thank Jens Stoye, Sven Rahmann, and Tobias Müller for fruitful discussions on
algorithmic and technical issues.
I appreciated very much the cleverly designed computing environment at the MPI for
Molecular Genetics maintained by Wilhelm Rüsing, Peter Marquardt, and others fromthe Computing Support group.
Furthermore I want to thank all friends, collaborators, and colleagues who gave support in any way, especially Antje Pleß, Christine Steinhoff, and all members of theAfter-Lunch-Walk Crew.
Financial support from the Deutsche Forschungsgemeinschaft and the Max-PlanckGesellschaft is gratefully acknowledged.
Finally, I want to thank the John-von-Neumann-Institute for Computing (NIC) inJülich and the Max-Planck-Institute for Molecular Genetics in Berlin for access to computing facilities.

111