Aleix Giménez Responsables: Blas Echebarría, Enrique Álvarez2013/201

Giménez, Aleix

English

UPCommons. Portal del coneixement obert de la UPC

Review of bioinformatics algorithms used in phylogenetic analysis and MSAAleix Giménez Grau∗ETSETB, Polytechnic University of Catalonia.(Dated: May 28, 2014)The problem of constructing phylogenetic trees is addressed from the point of view of bioinfor-matics, emphasizing their relation with multiple sequence alignments and presenting methods tosolve both problems. A C++ implementation of some algorithms was developed and applied to realnucleotide data downloaded from GenBank. By these means the phylogenetic tree of the Hominidaefamily is recovered with some minor errors.I. INTRODUCTIONBiology experimented a revolution in the middle of thepast century starting with the discovery of the structureof DNA in 1953 by Watson and Crick, and followed bygreat technical achievements, for instance the sequencingof the first complete genome in 1976. Nowadays exten-sive genetic information is freely available online fromGenBank1, and the main problem is how to treat datain order to obtain relevant results. Work has been donein this field in the last decades, resulting in the appear-ance of a new discipline, i.e. bioinformatics, that aims atusing computation in order to extract information frommolecular biological data.The study of evolution has been an important fieldsince the appearance of Darwin’s theory in 1859. Al-though comparative anatomy has been the main tool ofevolutionary biology during a long period bioinformat-ics has changed completely this methodology. Nowadaysspecies are no longer classified in terms of their anatomiccharacteristics but in terms of their genetic material,which has resulted in the concept of phylogenetic trees,i.e. evolutionary trees build in terms of genetic similarity.However, comparing the whole genome of species is stillnot possible, so shorter fragments must be used. Mito-chondrial DNA is widely used in phylogenetics because itis shorter and more stable than a complete genome, andit has been determined from many species, even extinctedones.In this paper, the methods for building phylogenetictrees based on bioinformatics will be reviewed. In orderto do so multiple sequence alignments (msa) must beanalyzed in detail and their link with phylogenetics mustbe stressed.The structure of the manuscript is the following: insection II the concept of pairwise alignment is definedand generalized to include multiple sequences. In sectionIII the main algorithms to find alignments and phyloge-netic trees are presented, i.e. Needleman and Wunschalgorithm for pairwise alignment, progressive alignmentfor msa, the Jukes-Cantor model for constructing dis-tance matrices and upgma and neighbor joining for phy-∗ aleixgim@gmail.comlogenetic tree building. Moreover, the implementation inC++ of these methods is summarized. In section IV theresults obtained applying the code to genetic data fromthe Hominidae family are shown and discussed. Finallyin section V conclusions are outlined briefly.II. ALIGNMENTS. GENERALITIESThe first problem that arises when dealing with geneticsequences of different species is how to compare them. Itis known that the origin of genetic diversity are muta-tions, that can take the form of substitution, deletion orinsertion of base pairs. An alignment is a way of arrang-ing the sequences of nucleotides or amino acids to identifyregions of similarity that may be a consequence of func-tional, structural, or evolutionary relationships betweenthe sequences2.We are going to discuss only nucleotide sequences, thatcan be thought as vectors of characters a, c, g, t. Giventwo sequences that can be of different lengths, an align-ment of them is an array which two rows are strings of thesame length with gaps in some positions, with the condi-tions that one column cannot contain two gaps and whenremoving the gaps the original strings must be recovered.Given the sequences aggtactatccagtca and atgctagattaat-gcaaggca a pairwise alignment that follows the previousconditions can be, for instanceag g tact at c ca g t caa t g ctag a t taatg caaggca(1)Given two sequences many alignments fulfill the previ-ous definition, so scoring functions are needed to definewhich is the best. A scoring function takes two argu-ments that can be characters representing nucleotides orgaps and it returns a real number. For DNA sequences itis customary to use the scoring function in table I: whenboth arguments are equal it is a match and the score isone; when the arguments are different it is a mismatchand the score is −µ; finally if one of the arguments isa gap it is an indel and the score is −σ. Although inthis work it is assumed the most simple situation withµ = σ = 1, there exists more advanced techniques, forinstance using distinct penalties for opening a gap or ex-tending it.2δ a c g t -a 1 -µ -µ -µ -σc -µ 1 -µ -µ -σg -µ -µ 1 -µ -σt -µ -µ -µ 1 -σ- -σ -σ -σ -σ ×TABLE I. Example of the simplest scoring function used innucleotide alignments. In this work it is assumed µ = σ = 1.When a scoring function is selected the total score ofthe alignment is defined as the sum of the score of allsites. Namely, if the alignment array is 2×l and its rowsare s′1j and s′2j the total score isSpw =l∑j=1δ(s′1j , s′2j). (2)A generalization of the previous score is needed if thealignment consists on N sequences arranged in the arrays′ij of the multiple sequence alignment. The sum of pairs(SP) score is defined as the sum of the score over allpairwise alignmentsSsp =N∑i=1N∑r=i+1 l∑j=1δ(s′ij , s′rj) . (3)The pairwise alignment problem consists on finding theoptimal position of the gaps in order to maximize thescore Spw, while the multiple sequence alignment prob-lem maximizes SSP . Although the formulation of bothproblems is very similar, the methods to solve them aresignificantly different, as will be seen in section III A.III. MATERIALS AND METHODSThe algorithms used in this work will be described inthis section, starting with techniques for pairwise andmultiple sequence alignment, and then presenting tree-building algorithms.A. Alignment techniquesOne approach to solve the pairwise alignment problemis the Needleman and Wunsch algorithm3, that reformu-lates it in terms of a directed acyclic graph (dag), andfinding its longest path is equivalent to solve the originalproblem. A dynamic programming algorithm is appliedto the dag so the alignment can be determined with com-putational complexity O(l2). This cost is not prohibitivefor most purposes with the current hardware of personalcomputers.In principle the multiple sequence alignment problemcould be solved with Needleman and Wunsch algorithm,but the computational complexity would be of O(2nln),which restrict the application of the exact algorithm toonly few short sequences. The problem has been shownto be NP-complete4, which roughly implies that an ef-ficient algorithm to solve it can not be found. This hasderived to the use of approximate or heuristic approachesto find msa, that make a trade-off between accuracy andefficiency.Progressive multiple alignment is one of the alogrithmsthat solves the previous problem approximately. Themain idea is that in the alignment of species with closephylogenetic relation few gaps need to be introduced,so these alignments are “better” than the ones renderedwith distant species. This suggests that a possible strat-egy to msa could be to pairwise align first the close se-quences and then progressively “align this alignments”starting from closer to more distant groups of species.Progressive alignment warranties that at least the align-ment of close species is carried properly, but convergenceto the optimal solution may not be reached.A more formal description of the algorithm is the fol-lowing: if the phylogenetic tree of the species is knownthe sequences can be progressively aligned following itstopology starting from leaves and arriving to the root.The phylogenetic tree previously used is called guidingtree, and can be obtained from a distance matrix method(see section III B). Most times a pairwise alignment oftwo arrays of multiple sequences must be done, and thisis accomplished using Needleman and Wunsch algorithmmaximizing SP-score of them. Using this approach thecomplexity of the algorithm is approximately O(2nl2)which is again a reasonable cost for most purposes.B. Distance-based tree-building methodsIn the previous section the relation between phy-logenetic trees and multiple sequence alignments wasglimpsed when a guide tree was used to find an approx-imate msa, but this link is tighter. A multiple align-ment somehow displays genetic information in such a waythat the evolution process that resulted in the studied se-quences can be determined easily. In fact, if two rows of amsa are almost identical they correspond to close speciesbut if there are a lot of mismatches then they correspondto distant species.When discussing the alignment of sequences scoringwas used as a measure of the similarity between them.However, in phylogenetics it is more useful the use ofdistances, being the simplest one the p-distance, whichis defined as the rate between the total number of mis-matches and indels and the length of the alignment:p =∑j p(s′1j , s′2j)l, p(a, b) ={1 if a 6= b0 if a = b(4)where a and b can be characters or gaps. The p-distanceis purely mathematical and do not make any biologicalassumption. More advanced distances use a model of evo-lution of DNA, being the simplest one the Jukes-Cantor3model, which assumes a constant mutation rate for allbase pairs and possibility of multiple mutations in onesite. The distance derived from it can be shown to be5d = −34ln(1− 4p3)(5)that for small distances (p ∼ 0) is equal to the p-distance,because it is unlikely that multiple mutations take placein the same site.Assume that multiple sequences of different species areknown, from what has been presented so far, species i andj could be pairwise aligned and the distance dij betweenthem could be computed with Jukes-Cantor model. Ifthis was done for all pairs, a matrix would be createdthat can be used in order to build phylogenetic trees,with the so called distance matrix methods.The simplest algorithm to do so is upgma6, that it-eratively clusters the closest two elements. When a newelement C is created from C1 and C2, the distance matrixis updated putting7d(C,C∗) = 0.5 d(C1, C∗) + 0.5 d(C2, C∗). (6)where C∗ is any other element. The principal drawbackof this approach is the fact that all branches have thesame length, so it is only valid under the assumption ofconstant evolution rate for all species.The distance method used in most situations isneighbor-joining8, that clusters elements following atrade-off between minimizing the distance between themand maximizing the distances with the other elements.Once a new element is created, formula (6) is used toupdate the distance matrix. It has been shown thatneighbor-joining produces trees that minimizes evolutionunder some particular circumstances9. This method pro-duces unrooted trees, that can be rooted by means of themid-point rule, which locates the root in the mid-point ofthe longest path of the tree, or with the outgroup methodthat assumes one of the species is the most distant tothe others and the root is positioned to obtain a treeconsistent with this hypothesis.C. ImplementationMost scientists use Python with extra modules forbioinformatics research; however a high level languageis not useful when the functions need to be programmedfrom scratch, because then they are significantly slower.With the use of a compiled language such as C++ this is nolonger a problem, but advanced features such as objectorientation can still be used. In order to derive the resultsof section IV a C++ code was developed and it has beenshared under a free license in Github10. It can be splitinto three main blocks: the first is responsible for generat-ing a pairwise alignment of two arrays of sequences usingthe Needleman and Wunsch algorithm and SP-scoringfunction, and can return the Jukes-Cantor distance be-tween two species; the second builds phylogenetic treesSpecie GB Code Specie GB CodeEuropean Human X90314 Chimp. schweinfurthii AF176722German neanderthal AF011222 Chimp. troglodytes AF176766Russian neanderthal AF254446 Chimp. vellerosus AF315498East. lowland gorilla AF050738 Chimp. verus AF176731East. mountain gorilla AF089820 Orangutan jari AF451964West. lowland gorilla AY079510 Orangutan puti AF451972TABLE II. Names of the species used to demonstrate the al-gorithms reviewed in the work. The Genbank code of theoriginal mitochondrial D-loop sequences is included for fur-ther reference; however, this sequences were modified in or-der to a assert that they corresponded to the same region ofDNA.with either upgma or neighbor-joining using a distancematrix that is automatically generated if the input isan array of all sequences; finally the third progressivelyaligns multiple sequences with the help of a guide treethat is created with the previous algorithms.IV. RESULTS AND DISCUSSIONThe main algorithms presented in this work were ap-plied to twelve sequences from species of the Hominidae(also known as great apes) family formed by humans,chimpanzees, gorillas and orangutans. A particular frag-ment of the mitochondrial D-loop was selected as thecomparing region, because it has been shown that D-loop nucleotides may be useful in studies of phylogeny invertebrates11. The considered species are listed in tableII together with their reference code in Genbank1.The genetic data was used to build a guiding tree fromwhich a multiple alignment was generated. This msa re-vealed that the sequences from table II do not representexactly the same segment of the D-loop, but some of themcontained extra nucleotides at the beginning or end. Thisincorrect regions could be easily removed from the orig-inal sequences comparing with the multiple alignment.The whole procedure was repeated with the correctedsequences and a new msa was generated which clearlycorresponded to an alignment of analogous sites of theDNA, as can be seen in figure 1.Once the sequences had been corrected the phyloge-netic trees were rebuild with upgma and neighbor-joiningmethods, exported to Newick format and printed usingthe Biopython package. The obtained results can be seenin figure 2, and it can be observed that the topology ofboth trees is the same, but different to the correct one, inwhich the branches of the orangutans and gorillas wouldbe swapped. It is worth noting that close species are clus-tered correctly and errors begin to appear when joiningdistant species, which is a very frequent phenomena inphylogenetics. As the topology of both trees is the samethe msa does not depend on the chosen guide tree, butit may not be the case in other situations.The correct topology could be obtained with neighbor-joining if the outgroup method was used assuming theorangutans were the most distant species. However, with4FIG. 1. First 160 nucleotides of the multiple sequence alignment obtained with a progressive method with SP-scoring and theupgma guiding tree of figure 2 (a).0.0 0.1 0.2 0.3 0.4branch length24681012taxa Gor_Western Gor_Eastern Gori_Rwanda Orangu_Jari Orangu_Puti Europ_human Neanthe_Rus Neanthe_Ger Chimp_Trogl Chimp_Schew Chimp_Verus Chimp_Velle0.0 0.1 0.2 0.3 0.4branch length24681012taxa Gor_Western Gori_Rwanda Gor_Eastern Orangu_Puti Orangu_Jari Europ_human Neanthe_Ger Neanthe_Rus Chimp_Velle Chimp_Verus Chimp_Schew Chimp_TroglFIG. 2. A matrix dij with all pairwise distances was generated with Needleman and Wunsch algorithm with Jukes-Cantordistance, and used to generate the trees (a) with upgma algorithm and (b) with neighbor-joining and mid-point rooting. Thegenetic data is the listed in table II but with the needed modifications to assure that the same sites of the mitochondrial D-loopwere compared.the used data there was not any reasonable argument tomake this assumption so this procedure was not carriedout. A fast look at the msa made it apparent that inthis specific region of the D-loop the most distant specieswere gorillas; this indicated that a more extensive regionon DNA should be probably used in order to recover thecorrect tree. Another possibility was that a better treebuilding method such as maximum parsimony, which hasnot been presented, could find the correct tree. Howeverthis hypothesis has not been tested in this work.V. CONCLUSIONSThis work has shown that genetic information is a pow-erful tool for phylogenetic studies, and that the adequatemethods can generate useful information that can be usedin multiple fields of biology. The bottleneck in the devel-opment of more advanced techniques is the lack of effi-cient algorithms to solve the problems, so a huge field ofresearch is open for mathematicians and computing sci-entists, who should find novel approaches that supersedethe current ones. However, the example at the end of thepaper has also shown that it is fundamental not to for-get the biological expectations and that it is importantto know that incorrect results may be obtained so theyneed to be contrasted properly.[1] D. A. Benson, I. Karsch-Mizrachi, D. J. Lipman, J. Os-tell, and D. L. Wheeler, Nucleic Acids Research 33, D34(2005).[2] D. W. Mount, Bioinformatics: Sequence and GenomeAnalysis (Cold Spring Harbor Laboratory Press, 2004).[3] S. B. Needleman and C. Wunsch, Journal of MolecularBiology 48, 443 (1970).[4] L. Wang and T. Jiang, Journal of computational biology51, 337 (1994).[5] T. H. Jukes and C. R. Cantor, Evolution of ProteinMolecules, edited by H. N. Munro (Academy Press,1969).[6] R. R. Sokal and C. D. Michener, University of KansasScientific Bulletin 28, 1409 (1958).[7] H. Böckenhauer and D. Bongartz, Algorithmic Aspectsof Bioinformatics, Natural Computing Series (Springer,2007).[8] N. Saitou and M. Nei, Molecular Biology and Evolution4, 406 (1987).[9] O. Gascuel and M. Steel, Molecular Biology and Evolu-tion 23, 1997 (2006).[10] https://github.com/aleixgg/Bioinformatics CPP.[11] A. Larizza, G. Pesole, A. Reyes, E. Sbis, and C. Saccone,Journal of Molecular Evolution 54, 145 (2002).

Review of bioinformatics algorithms used in phylogenetic analysis and MSA

https://upcommons.upc.edu/bitstream/2117/190534/1/bioinformaticaaleixgg1.pdf

Review of bioinformatics algorithms used in phylogenetic analysis and MSA

Abstract

Similar works

Full text

Available Versions

UPCommons. Portal del coneixement obert de la UPC