-
Notifications
You must be signed in to change notification settings - Fork 0
/
emma.html
1719 lines (1363 loc) · 71.2 KB
/
emma.html
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
<HTML>
<HEAD>
<TITLE>
EMBOSS: emma
</TITLE>
</HEAD>
<BODY BGCOLOR="#FFFFFF" text="#000000">
<table align=center border=0 cellspacing=0 cellpadding=0>
<tr><td valign=top>
<A HREF="/" ONMOUSEOVER="self.status='Go to the EMBOSS home page';return true"><img border=0 src="emboss_icon.jpg" alt="" width=150 height=48></a>
</td>
<td align=left valign=middle>
<b><font size="+6">
emma
</font></b>
</td></tr>
</table>
<br>
<p>
<H2>
Function
</H2>
Multiple alignment program - interface to ClustalW program
<H2>
Description
</H2>
EMMA calculates the multiple alignment of nucleic acid or protein
sequences according to the method of Thompson, J.D., Higgins, D.G.
and Gibson, T.J. (1994).
<p>
This is an interface to the ClustalW distribution.
<h3>The basic alignment method</h3>
The basic multiple alignment algorithm consists of three main stages: 1) all
pairs of sequences are aligned separately in order to calculate a distance matrix
giving the divergence of each pair of sequences; 2) a guide tree is calculated
from the distance matrix; 3) the sequences are progressively aligned according
to the branching order in the guide tree. An example using 7 globin
sequences of known tertiary structure (25) is given in figure 1.
<h4>1) The distance matrix/pairwise alignments</h4>
In the original CLUSTAL programs, the pairwise distances were calculated
using a fast approximate method (22). This allows very large numbers of
sequences to be aligned, even on a microcomputer. The scores are calculated
as the number of k-tuple matches (runs of identical residues, typically 1 or 2
long for proteins or 2 to 4 long for nucleotide sequences) in the best alignment
between two sequences minus a fixed penalty for every gap. We now offer a
choice between this method and the slower but more accurate scores from full
dynamic programming alignments using two gap penalties (for opening or
extending gaps) and a full amino acid weight matrix. These scores are
calculated as the number of identities in the best alignment divided by the
number of residues compared (gap positions are excluded). Both of these
scores are initially calculated as percent identity scores and are converted to
distances by dividing by 100 and subtracting from 1.0 to give number of
differences per site. We do not correct for multiple substitutions in these
initial distances. In figure 1 we give the 7x7 distance matrix between the 7
globin sequences calculated using the full dynamic programming method.
<h4>2) The guide tree</h4>
The trees used to guide the final multiple alignment process are calculated
from the distance matrix of step 1 using the Neighbour-Joining method (21).
This produces unrooted trees with branch lengths proportional to estimated
divergence along each branch. The root is placed by a "mid-point" method
(15) at a position where the means of the branch lengths on either side of the
root are equal. These trees are also used to derive a weight for each sequence
(15). The weights are dependent upon the distance from the root of the tree
but sequences which have a common branch with other sequences share the
weight derived from the shared branch. In the example in figure 1, the
leghaemoglobin (Lgb2_Luplu) gets a weight of 0.442 which is equal to the
length of the branch from the root to it. The Human beta globin
(Hbb_Human) gets a weight consisting of the length of the branch leading to
it that is not shared with any other sequences (0.081) plus half the length of
the branch shared with the horse beta globin (0.226/2) plus one quarter the
length of the branch shared by all four haemoglobins (0.061/4) plus one fifth
the branch shared between the haemoglobins and the myoglobin (0.015/5)
plus one sixth the branch leading to all the vertebrate globins (0.062). This
sums to a total of 0.221. By contrast, in the normal progressive alignment
algorithm, all sequences would be equally weighted. The rooted tree with
branch lengths and sequence weights for the 7 globins is given in figure 1.
<h4>3) Progressive alignment</h4>
The basic procedure at this stage is to use a series of pairwise alignments to
align larger and larger groups of sequences, following the branching order in
the guide tree. You proceed from the tips of the rooted tree towards the root.
<p>
In the globin example in figure 1 you align the sequences in the following
order: human vs. horse beta globin; human vs. horse alpha globin; the 2
alpha globins vs. the 2 beta globins; the myoglobin vs. the haemoglobins; the
cyanohaemoglobin vs the haemoglobins plus myoglobin; the leghaemoglobin
vs. all the rest. At each stage a full dynamic programming (26,27) algorithm is
used with a residue weight matrix and penalties for opening and extending
gaps. Each step consists of aligning two existing alignments or sequences.
Gaps that are present in older alignments remain fixed. In the basic
algorithm, new gaps that are introduced at each stage get full gap opening and
extension penalties, even if they are introduced inside old gap positions (see
the section on gap penalties below for modifications to this rule). In order to
calculate the score between a position from one sequence or alignment and
one from another, the average of all the pairwise weight matrix scores from
the amino acids in the two sets of sequences is used i.e. if you align 2
alignments with 2 and 4 sequences respectively, the score at each position is
the average of 8 (2x4) comparisons. This is illustrated in figure 2. If either set
of sequences contains one or more gaps in one of the positions being
considered, each gap versus a residue is scored as zero. The default amino
acid weight matrices we use are rescored to have only positive values.
Therefore, this treatment of gaps treats the score of a residue versus a gap as
having the worst possible score. When sequences are weighted (see
improvements to progressive alignment, below), each weight matrix value is
multiplied by the weights from the 2 sequences, as illustrated in figure 2.
<h4>Improvements to progressive alignment</h4>
All of the remaining modifications apply only to the final progressive
alignment stage. Sequence weighting is relatively straightforward and is
already widely used in profile searches (15,16). The treatment of gap penalties
is more complicated. Initial gap penalties are calculated depending on the
weight matrix, the similarity of the sequences, and the length of the
sequences. Then, an attempt is made to derive sensible local gap opening
penalties at every position in each pre-aligned group of sequences that will
vary as new sequences are added. The use of different weight matrices as the
alignment progresses is novel and largely by-passes the problem of initial
choice of weight matrix. The final modification allows us to delay the
addition of very divergent sequences until the end of the alignment process
when all of the more closely related sequences have already been aligned.
<h3>Sequence weighting</h3>
Sequence weights are calculated directly from the guide tree. The weights
are normalised such that the biggest one is set to 1.0 and the rest are all less
than one. Groups of closely related sequences receive lowered weights
because they contain much duplicated information. Highly divergent
sequences without any close relatives receive high weights. These weights
are used as simple multiplication factors for scoring positions from different
sequences or prealigned groups of sequences. The method is illustrated in
figure 2. In the globin example in figure 1, the two alpha globins get
downweighted because they are almost duplicate sequences (as do the two
beta globins); they receive a combined weight of only slightly more than if a
single alpha globin was used.
<h3>Initial gap penalties</h3>
Initially, two gap penalties are used: a gap opening penalty (GOP) which gives
the cost of opening a new gap of any length and a gap extension penalty (GEP)
which gives the cost of every item in a gap. Initial values can be set by the
user from a menu. The software then automatically attempts to choose
appropriate gap penalties for each sequence alignment, depending on the
following factors.
<h4>1) Dependence on the weight matrix</h4>
It has been shown (16,28) that varying the gap penalties used with different
weight matrices can improve the accuracy of sequence alignments. Here, we
use the average score for two mismatched residues (ie. off-diagonal values in
the matrix) as a scaling factor for the GOP.
<h4>2) Dependence on the similarity of the sequences</h4>
The percent identity of the two (groups of) sequences to be aligned is used to
increase the GOP for closely related sequences and decrease it for more
divergent sequences on a linear scale.
<h4>3) Dependence on the lengths of the sequences</h4>
The scores for both true and false sequence alignments grow with the length
of the sequences. We use the logarithm of the length of the shorter sequence
to increase the GOP with sequence length.
<p>
Using these three modifications, the initial GOP calculated by the program is:
<p>
GOP->(GOP+log(MIN(N,M))) * (average residue mismatch score) * (percent identity scaling factor)
<br>
where N, M are the lengths of the two sequences.
<p>
<h4>4) Dependence on the difference in the lengths of the sequences</h4>
The GEP is modified depending on the difference between the lengths of the
two sequences to be aligned. If one sequence is much shorter than the other,
the GEP is increased to inhibit too many long gaps in the shorter sequence.
The initial GEP calculated by the program is:
<p>
GEP -> GEP*(1.0+|log(N/M)|)
<br>
where N, M are the lengths of the two sequences.
<h3>Position-specific gap penalties</h3>
In most dynamic programming applications, the initial gap opening and
extension penalties are applied equally at every position in the sequence,
regardless of the location of a gap, except for terminal gaps which are usually
allowed at no cost. In CLUSTAL W, before any pair of sequences or
prealigned groups of sequences are aligned, we generate a table of gap opening
penalties for every position in the two (sets of) sequences. An example is
shown in figure 3. We manipulate the initial gap opening penalty in a
position specific manner, in order to make gaps more or less likely at different
positions.
<p>
The local gap penalty modification rules are applied in a hierarchical manner.
<p>
The exact details of each rule are given below. Firstly, if there is a gap at a
position, the gap opening and gap extension penalties are lowered; the other
rules do not apply. This makes gaps more likely at positions where there are
already gaps. If there is no gap at a position, then the gap opening penalty is
increased if the position is within 8 residues of an existing gap. This
discourages gaps that are too close together. Finally, at any position within a
run of hydrophilic residues, the penalty is decreased. These runs usually
indicate loop regions in protein structures. If there is no run of hydrophilic
residues, the penalty is modified using a table of residue specific gap
propensities (12). These propensities were derived by counting the frequency
of each residue at either end of gaps in alignments of proteins of known
structure. An illustration of the application of these rules from one part of
the globin example, in figure 1, is given in figure 3.
<h4>1) Lowered gap penalties at existing gaps</h4>
If there are already gaps at a position, then the GOP is reduced in proportion
to the number of sequences with a gap at this position and the GEP is lowered
by a half. The new gap opening penalty is calculated as:
<p>
GOP -> GOP*0.3*(no. of sequences without a gap/no. of sequences).
<h4>2) Increased gap penalties near existing gaps</h4>
If a position does not have any gaps but is within 8 residues of an existing gap,
the GOP is increased by:
<p>
GOP -> GOP*(2+((8-distance from gap)*2)/8)
<h4>3) Reduced gap penalties in hydrophilic stretches</h4>
Any run of 5 hydrophilic residues is considered to be a hydrophilic stretch.
The residues that are to be considered hydrophilic may be set by the user but
are conservatively set to D, E, G, K, N, Q, P, R or S by default. If, at any
position, there are no gaps and any of the sequences has such a stretch, the
GOP is reduced by one third.
<h4>4) Residue specific penalties</h4>
If there is no hydrophilic stretch and the position does not contain any gaps,
then the GOP is multiplied by one of the 20 numbers in table 1, depending on
the residue. If there is a mixture of residues at a position, the multiplication
factor is the average of all the contributions from each sequence.
<h3>Weight matrices</h3>
Two main series of weight matrices are offered to the user: the Dayhoff PAM
series (3) and the BLOSUM series (4). The default is the BLOSUM series. In
each case, there is a choice of matrix ranging from strict ones, useful for
comparing very closely related sequences to very "soft" ones that are useful
for comparing very distantly related sequences. Depending on the distance
between the two sequences or groups of sequences to be compared, we switch
between 4 different matrices. The distances are measured directly from the
guide tree. The ranges of distances and tables used with the PAM series of
matrices is: 80-100%:PAM20, 60-80%:PAM60, 40-60%:PAM120, 0-40%:PAM350.
The range used with the BLOSUM series is:80-100%:BLOSUM80,
60-80%:BLOSUM62, 30-60%:BLOSUM45, 0-30%:BLOSUM30.
<h3>Divergent sequences</h3>
The most divergent sequences (most different, on average from all of the
other sequences) are usually the most difficult to align correctly. It is
sometimes better to delay the incorporation of these sequences until all of the
more easily aligned sequences are merged first. This may give a better chance
of correctly placing the gaps and matching weakly conserved positions against
the rest of the sequences. A choice is offered to set a cut off (default is 40%
identity or less with any other sequence) that will delay the alignment of the
divergent sequences until all of the rest have been aligned.
<h2>Software and Algorithms</h2>
<h3>Dynamic Programming</h3>
The most demanding part of the multiple alignment strategy, in terms of
computer processing and memory usage, is the alignment of two (groups of)
sequences at each step in the final progressive alignment. To make it
possible to align very long sequences (e.g. dynein heavy chains at ~ 5,000
residues) in a reasonable amount of memory, we use the memory efficient
dynamic programming algorithm of Myers and Miller (26). This sacrifices
some processing time but makes very large alignments practical in very little
memory. One disadvantage of this algorithm is that it does not allow
different gap opening and extension penalties at each position. We have
modified the algorithm so as to allow this and the details are described in a
separate paper (27).
<h3>Alignment to an alignment</h3>
Profile alignment is used to align two existing alignments (either of which
may consist of just one sequence) or to add a series of new sequences to an
existing alignment. This is useful because one may wish to build up a
multiple alignment gradually, choosing different parameters manually, or
correcting intermediate errors as the alignment proceeds. Often, just a few
sequences cause misalignments in the progressive algorithm and these can be
removed from the process and then added at the end by profile alignment. A
second use is where one has a high quality reference alignment and wishes to
keep it fixed while adding new sequences automatically.
<!-- clustalw.doc documentation -->
<h3>Terminal Gaps</h3>
In the original Clustal V program, terminal gaps were penalised the same
as all other gaps. This caused some ugly side effects e.g.
<p>
<pre>
acgtacgtacgtacgt acgtacgtacgtacgt
a----cgtacgtacgt gets the same score as ----acgtacgtacgt
</pre>
<p>
NOW, terminal gaps are free. This is better on average and stops silly
effects like single residues jumping to the edge of the alignment. However,
it is not perfect. It does mean that if there should be a gap near the end
of the alignment, the program may be reluctant to insert it i.e.
<p>
<pre>
cccccgggccccc cccccgggccccc
ccccc---ccccc may be considered worse (lower score) than cccccccccc---
</pre>
<p>
In the right hand case above, the terminal gap is free and may score higher
than the laft hand alignment. This can be prevented by lowering the gap
opening and extension penalties. It is difficult to get this right all the
time. Please watch the ends of your alignments.
<h3>Speed of the initial (pairwise) alignments (fast approximate/slow accurate)</h3>
By default, the initial pairwise alignments are now carried out using a full
dynamic programming algorithm. This is more accurate than the older hash/
k-tuple based alignments (Wilbur and Lipman) but is MUCH slower. On a fast
workstation you may not notice but on a slow box, the difference is extreme.
You can set the alignment method from the menus easily to the older, faster
method.
<h3>Delaying alignment of distant sequences</h3>
The user can set a cut off to delay the alignment of the most divergent
sequences in a data set until all other sequences have been aligned. By
default, this is set to 40% which means that if a sequence is less than 40%
identical to any other sequence, its alignment will be delayed.
<h3>Iterative realignment/Reset gaps between alignments</h3>
By default, if you align a set of sequences a second time (e.g. with changed
gap penalties), the gaps from the first alignment are discarded. You can
set this from the menus so that older gaps will be kept between alignments,
This can sometimes give better alignments by keeping the gaps (do not reset
them) and doing the full multiple alignment a second time. Sometimes, the
alignment will converge on a better solution; sometimes the new alignment will
be the same as the first. There can be a strange side effect: you can get
columns of nothing but gaps introduced.
<p>
Any gaps that are read in from the input file are always kept, regardless
of the setting of this switch. If you read in a full multiple alignment, the "reset
gaps" switch has no effect. The old gaps will remain and if you carry out
a multiple alignment, any new gaps will be added in. If you wish to carry out
a full new alignment of a set of sequences that are already aligned in a file
you must input the sequences without gaps.
<h3>Profile alignment</h3>
By profile alignment, we simply mean the alignment of old alignments/sequences.
In this context, a profile is just an existing alignment (or even a set of
unaligned sequences; see below). This allows you to
read in an old alignment (in any of the allowed input formats) and align
one or more new sequences to it. From the profile alignment menu, you
are allowed to read in 2 profiles. Either profile can be a full alignment
OR a single sequence. In the simplest mode, you simply align the two profiles
to each other. This is useful if you want to gradually build up a full
multiple alignment.
<p>
A second option is to align the sequences from the second profile, one at
a time to the first profile. This is done, taking the underlying tree between
the sequences into account. This is useful if you have a set of new sequences
(not aligned) and you wish to add them all to an older alignment.
<h2>Changes to the phylogentic tree calculations and some hints</h2>
<h3>Improved distance calculations for protein trees</h3>
The phylogenetic trees in Clustal W (the real trees that you calculate
AFTER alignment; not the guide trees used to decide the branching order
for multiple alignment) use the Neighbor-Joining method of Saitou and
Nei based on a matrix of "distances" between all sequences. These distances
can be corrected for "multiple hits". This is normal practice when accurate
trees are needed. This correction stretches distances (especially large ones)
to try to correct for the fact that OBSERVED distances (mean number of
differences per site) greatly underestimate the actual number that happened
during evolution.
<p>
In Clustal V we used a simple formula to convert an observed distance to one
that is corrected for multiple hits. The observed distance is the mean number
of differences per site in an alignment (ignoring sites with a gap) and is
therefore always between 0.0 (for ientical sequences) an 1.0 (no residues the
same at any site). These distances can be multiplied by 100 to give percent
difference values. 100 minus percent difference gives percent identity.
The formula we use to correct for multiple hits is from Motoo Kimura
(Kimura, M. The neutral Theory of Molecular Evolution, Camb.Univ.Press, 1983,
page 75) and is:
<p>
K = -Ln(1 - D - (D.D)/5)
<br>
where D is the observed distance and K is
corrected distance.
<p>
This formula gives mean number of estimated substitutions per site and, in
contrast to D (the observed number), can be greater than 1 i.e. more than
one substitution per site, on average. For example, if you observe 0.8
differences per site (80% difference; 20% identity), then the above formula
predicts that there have been 2.5 substitutions per site over the course
of evolution since the 2 sequences diverged. This can also be expressed in
PAM units by multiplying by 100 (mean number of substitutions per 100 residues).
The PAM scale of evolution and its derivation/calculation comes from the
work of Margaret Dayhoff and co workers (the famous Dayhoff PAM series
of weight matrices also came from this work). Dayhoff et al constructed
an elaborate model of protein evolution based on observed frequencies
of substitution between very closely related proteins. Using this model,
they derived a table relating observed distances to predicted PAM distances.
Kimura's formula, above, is just a "curve fitting" approximation to this table.
It is very accurate in the range 0.75 > D > 0.0 but becomes increasingly
unaccurate at high D (>0.75) and fails completely at around D = 0.85.
<p>
To circumvent this problem, we calculated all the values for K corresponding
to D above 0.75 directly using the Dayhoff model and store these in an
internal table, used by Clustal W. This table is declared in the file dayhoff.h and
gives values of K for all D between 0.75 and 0.93 in intervals of 0.001 i.e.
for D = 0.750, 0.751, 0.752 ...... 0.929, 0.930. For any observed D
higher than 0.930, we arbitrarily set K to 10.0. This sounds drastic but
with real sequences, distances of 0.93 (less than 7% identity) are rare.
If your data set includes sequences with this degree of divergence, you
will have great difficulty getting accurate trees by ANY method; the alignment
itself will be very difficult (to construct and to evaluate).
<p>
There are some important
things to note. Firstly, this formula works well if your sequences are
of average amino acid composition and if the amino acids substitute according
to the original Dayhoff model. In other cases, it may be misleading. Secondly,
it is based only on observed percent distance i.e. it does not DIRECTLY
take conservative substitutions into account. Thirdly, the error on the
estimated PAM distances may be VERY great for high distances; at very high
distance (e.g. over 85%) it may give largely arbitrary corrected distances.
In most cases, however, the correction is still worth using; the trees will
be more accurate and the branch lengths will be more realistic.
<p>
A far more sophisticated distance correction based on a full Dayhoff
model which DOES take conservative substitutions and actual amino acid
composition into account, may be found in the PROTDIST program of the
PHYLIP package. For serious tree makers, this program is highly recommended.
<h3>TWO NOTES ON BOOTSTRAPPING...</h3>
When you use the BOOTSTRAP in Clustal W to estimate the reliability of parts
of a tree, many of the uncorrected distances may randomly exceed the arbitrary cut
off of 0.93 (sequences only 7% identical) if the sequences are distantly
related. This will happen randomly i.e. even if none of the pairs of
sequences are less than 7% identical, the bootstrap samples may contain pairs
of sequences that do exceed this cut off.
If this happens, you will be warned. In practice, this can
happen with many data sets. It is not a serious problem if it happens rarely.
If it does happen (you are warned when it happens and told how often the
problem occurs), you should consider removing the most distantly
related sequences and/or using the PHYLIP package instead.
<p>
A further problem arises in almost exactly the opposite situation: when
you bootstrap a data set which contains 3 or more sequences that are identical
or almost identical. Here, the sets of identical sequences should be shown
as a multifurcation (several sequences joing at the same part of the tree).
Because the Neighbor-Joining method only gives strictly dichotomous trees
(never more than 2 sequences join at one time), this cannot be exactly
represented. In practice, this is NOT a problem as there will be some
internal branches of zero length seperating the sequences. If you
display the tree with all branch lengths, you will still see a multifurcation.
However, when you bootstrap
the tree, only the branching orders are stored and counted. In the case
of multifurcations, the exact branching order is arbitrary but the program
will always get the same branching order, depending only on the input order
of the sequences. In practice, this is only a problem in situations where
you have a set of sequences where all of them are VERY similar. In this case,
you can find very high support for some groupings which will disappear if you
run the analysis with a different input order. Again, the PHYLIP package
deals with this by offering a JUMBLE option to shuffle the input order
of your sequences between each bootstrap sample.
<H2>
Usage
</H2>
<b>Here is a sample session with emma</b>
<p>
<p>
<table width="90%"><tr><td bgcolor="#CCFFFF"><pre>
% <b>emma </b>
Multiple alignment program - interface to ClustalW program
Input (gapped) sequence(s): <b>globins.fasta</b>
output sequence set [hbb_human.aln]: <b></b>
Dendrogram (tree file) from clustalw output file [hbb_human.dnd]: <b></b>
CLUSTAL W (1.83) Multiple Sequence Alignments
Sequence type explicitly set to Protein
Sequence format is Pearson
Sequence 1: HBB_HUMAN 146 aa
Sequence 2: HBB_HORSE 146 aa
Sequence 3: HBA_HUMAN 141 aa
Sequence 4: HBA_HORSE 141 aa
Sequence 5: MYG_PHYCA 153 aa
Sequence 6: GLB5_PETMA 149 aa
Sequence 7: LGB2_LUPLU 153 aa
Start of Pairwise alignments
Aligning...
Sequences (1:2) Aligned. Score: 83
Sequences (1:3) Aligned. Score: 43
Sequences (1:4) Aligned. Score: 42
Sequences (1:5) Aligned. Score: 24
Sequences (1:6) Aligned. Score: 21
Sequences (1:7) Aligned. Score: 14
Sequences (2:3) Aligned. Score: 41
Sequences (2:4) Aligned. Score: 43
Sequences (2:5) Aligned. Score: 24
Sequences (2:6) Aligned. Score: 19
Sequences (2:7) Aligned. Score: 15
Sequences (3:4) Aligned. Score: 87
Sequences (3:5) Aligned. Score: 26
Sequences (3:6) Aligned. Score: 29
Sequences (3:7) Aligned. Score: 16
Sequences (4:5) Aligned. Score: 26
Sequences (4:6) Aligned. Score: 27
Sequences (4:7) Aligned. Score: 12
Sequences (5:6) Aligned. Score: 21
Sequences (5:7) Aligned. Score: 7
Sequences (6:7) Aligned. Score: 11
Guide tree file created: [12345678A]
Start of Multiple Alignment
There are 6 groups
Aligning...
Group 1: Sequences: 2 Score:2194
Group 2: Sequences: 2 Score:2165
Group 3: Sequences: 4 Score:960
Group 4: Delayed
Group 5: Delayed
Group 6: Delayed
Sequence:5 Score:865
Sequence:6 Score:797
Sequence:7 Score:1044
Alignment Score 4164
GCG-Alignment file created [12345678A]
</pre></td></tr></table><p>
<p>
<a href="#input.1">Go to the input files for this example</a><br><a href="#output.1">Go to the output files for this example</a><p><p>
<H2>
Command line arguments
</H2>
<table CELLSPACING=0 CELLPADDING=3 BGCOLOR="#f5f5ff" ><tr><td>
<pre>
Standard (Mandatory) qualifiers:
[-sequence] seqall (Gapped) sequence(s) filename and optional
format, or reference (input USA)
[-outseq] seqoutset [<sequence>.<format>] Sequence set filename
and optional format (output USA)
[-dendoutfile] outfile [*.emma] Dendrogram (tree file) from
clustalw output file
Additional (Optional) qualifiers (* if not always prompted):
-onlydend toggle [N] Only produce dendrogram file
* -dend toggle [N] Do alignment using an old dendrogram
* -dendfile infile Dendrogram (tree file) from clustalw file
(optional)
* -pwmatrix menu [b] The scoring table which describes the
similarity of each amino acid to each other.
There are three 'in-built' series of weight
matrices offered. Each consists of several
matrices which work differently at different
evolutionary distances. To see the exact
details, read the documentation. Crudely, we
store several matrices in memory, spanning
the full range of amino acid distance (from
almost identical sequences to highly
divergent ones). For very similar sequences,
it is best to use a strict weight matrix
which only gives a high score to identities
and the most favoured conservative
substitutions. For more divergent sequences,
it is appropriate to use 'softer' matrices
which give a high score to many other
frequent substitutions.
1) BLOSUM (Henikoff). These matrices appear
to be the best available for carrying out
data base similarity (homology searches).
The matrices used are: Blosum80, 62, 45 and
30.
2) PAM (Dayhoff). These have been extremely
widely used since the late '70s. We use the
PAM 120, 160, 250 and 350 matrices.
3) GONNET . These matrices were derived
using almost the same procedure as the
Dayhoff one (above) but are much more up to
date and are based on a far larger data set.
They appear to be more sensitive than the
Dayhoff series. We use the GONNET 40, 80,
120, 160, 250 and 350 matrices.
We also supply an identity matrix which
gives a score of 1.0 to two identical amino
acids and a score of zero otherwise. This
matrix is not very useful. (Values: b
(blosum); p (pam); g (gonnet); i (id); o
(own))
* -pwdnamatrix menu [i] The scoring table which describes the
scores assigned to matches and mismatches
(including IUB ambiguity codes). (Values: i
(iub); c (clustalw); o (own))
* -pairwisedatafile infile Comparison matrix file (optional)
* -matrix menu [b] This gives a menu where you are offered
a choice of weight matrices. The default for
proteins is the PAM series derived by
Gonnet and colleagues. Note, a series is
used! The actual matrix that is used depends
on how similar the sequences to be aligned
at this alignment step are. Different
matrices work differently at each
evolutionary distance.
There are three 'in-built' series of weight
matrices offered. Each consists of several
matrices which work differently at different
evolutionary distances. To see the exact
details, read the documentation. Crudely, we
store several matrices in memory, spanning
the full range of amino acid distance (from
almost identical sequences to highly
divergent ones). For very similar sequences,
it is best to use a strict weight matrix
which only gives a high score to identities
and the most favoured conservative
substitutions. For more divergent sequences,
it is appropriate to use 'softer' matrices
which give a high score to many other
frequent substitutions.
1) BLOSUM (Henikoff). These matrices appear
to be the best available for carrying out
data base similarity (homology searches).
The matrices used are: Blosum80, 62, 45 and
30.
2) PAM (Dayhoff). These have been extremely
widely used since the late '70s. We use the
PAM 120, 160, 250 and 350 matrices.
3) GONNET . These matrices were derived
using almost the same procedure as the
Dayhoff one (above) but are much more up to
date and are based on a far larger data set.
They appear to be more sensitive than the
Dayhoff series. We use the GONNET 40, 80,
120, 160, 250 and 350 matrices.
We also supply an identity matrix which
gives a score of 1.0 to two identical amino
acids and a score of zero otherwise. This
matrix is not very useful. Alternatively,
you can read in your own (just one matrix,
not a series). (Values: b (blosum); p (pam);
g (gonnet); i (id); o (own))
* -dnamatrix menu [i] This gives a menu where a single matrix
(not a series) can be selected. (Values: i
(iub); c (clustalw); o (own))
* -mamatrixfile infile Comparison matrix file (optional)'
-[no]slow toggle [Y] A distance is calculated between every
pair of sequences and these are used to
construct the dendrogram which guides the
final multiple alignment. The scores are
calculated from separate pairwise
alignments. These can be calculated using 2
methods: dynamic programming (slow but
accurate) or by the method of Wilbur and
Lipman (extremely fast but approximate).
The slow-accurate method is fine for short
sequences but will be VERY SLOW for many
(e.g. >100) long (e.g. >1000 residue)
sequences.
* -pwgapopen float [10.0] The penalty for opening a gap in the
pairwise alignments. (Number 0.000 or more)
* -pwgapextend float [0.1] The penalty for extending a gap by 1
residue in the pairwise alignments. (Number
0.000 or more)
* -ktup integer [1 for protein, 2 for nucleic] This is the
size of exactly matching fragment that is
used. INCREASE for speed (max= 2 for
proteins; 4 for DNA), DECREASE for
sensitivity. For longer sequences (e.g.
>1000 residues) you may need to increase the
default. (integer from 0 to 4)
* -gapw integer [3 for protein, 5 for nucleic] This is a
penalty for each gap in the fast alignments.
It has little affect on the speed or
sensitivity except for extreme values.
(Positive integer)
* -topdiags integer [5 for protein, 4 for nucleic] The number of
k-tuple matches on each diagonal (in an
imaginary dot-matrix plot) is calculated.
Only the best ones (with most matches) are
used in the alignment. This parameter
specifies how many. Decrease for speed;
increase for sensitivity. (Positive integer)
* -window integer [5 for protein, 4 for nucleic] This is the
number of diagonals around each of the
'best' diagonals that will be used. Decrease
for speed; increase for sensitivity.
(Positive integer)
* -nopercent boolean [N] Fast pairwise alignment: similarity
scores: suppresses percentage score
-gapopen float [10.0] The penalty for opening a gap in the
alignment. Increasing the gap opening
penalty will make gaps less frequent.
(Positive foating point number)
-gapextend float [5.0] The penalty for extending a gap by 1
residue. Increasing the gap extension
penalty will make gaps shorter. Terminal
gaps are not penalised. (Positive foating
point number)
-[no]endgaps boolean [Y] End gap separation: treats end gaps just
like internal gaps for the purposes of
avoiding gaps that are too close (set by
'gap separation distance'). If you turn this
off, end gaps will be ignored for this
purpose. This is useful when you wish to
align fragments where the end gaps are not
biologically meaningful.
-gapdist integer [8] Gap separation distance: tries to
decrease the chances of gaps being too close
to each other. Gaps that are less than this
distance apart are penalised more than
other gaps. This does not prevent close
gaps; it makes them less frequent, promoting
a block-like appearance of the alignment.
(Positive integer)
* -norgap boolean [N] Residue specific penalties: amino acid
specific gap penalties that reduce or
increase the gap opening penalties at each
position in the alignment or sequence. As an
example, positions that are rich in glycine
are more likely to have an adjacent gap
than positions that are rich in valine.
* -hgapres string [GPSNDQEKR] This is a set of the residues
'considered' to be hydrophilic. It is used
when introducing Hydrophilic gap penalties.
(Any string is accepted)
* -nohgap boolean [N] Hydrophilic gap penalties: used to
increase the chances of a gap within a run
(5 or more residues) of hydrophilic amino
acids; these are likely to be loop or random
coil regions where gaps are more common.
The residues that are 'considered' to be
hydrophilic are set by '-hgapres'.
-maxdiv integer [30] This switch, delays the alignment of
the most distantly related sequences until
after the most closely related sequences
have been aligned. The setting shows the
percent identity level required to delay the
addition of a sequence; sequences that are
less identical than this level to any other
sequences will be aligned later. (Integer
from 0 to 100)
Advanced (Unprompted) qualifiers: (none)
Associated qualifiers:
"-sequence" associated qualifiers
-sbegin1 integer Start of each sequence to be used
-send1 integer End of each sequence to be used
-sreverse1 boolean Reverse (if DNA)
-sask1 boolean Ask for begin/end/reverse
-snucleotide1 boolean Sequence is nucleotide
-sprotein1 boolean Sequence is protein
-slower1 boolean Make lower case
-supper1 boolean Make upper case
-sformat1 string Input sequence format
-sdbname1 string Database name
-sid1 string Entryname
-ufo1 string UFO features
-fformat1 string Features format
-fopenfile1 string Features file name
"-outseq" associated qualifiers
-osformat2 string Output seq format
-osextension2 string File name extension
-osname2 string Base file name
-osdirectory2 string Output directory
-osdbname2 string Database name to add
-ossingle2 boolean Separate file for each entry
-oufo2 string UFO features
-offormat2 string Features format
-ofname2 string Features file name
-ofdirectory2 string Output directory
"-dendoutfile" associated qualifiers
-odirectory3 string Output directory
General qualifiers:
-auto boolean Turn off prompts
-stdout boolean Write standard output
-filter boolean Read standard input, write standard output
-options boolean Prompt for standard and additional values
-debug boolean Write debug output to program.dbg
-verbose boolean Report some/full command line options
-help boolean Report command line options. More
information on associated and general
qualifiers can be found with -help -verbose
-warning boolean Report warnings
-error boolean Report errors
-fatal boolean Report fatal errors
-die boolean Report dying program messages
</pre>
</td></tr></table>
<P>
<table border cellspacing=0 cellpadding=3 bgcolor="#ccccff">
<tr bgcolor="#FFFFCC">
<th align="left" colspan=2>Standard (Mandatory) qualifiers</th>
<th align="left">Allowed values</th>
<th align="left">Default</th>
</tr>
<tr>
<td>[-sequence]<br>(Parameter 1)</td>
<td>(Gapped) sequence(s) filename and optional format, or reference (input USA)</td>
<td>Readable sequence(s)</td>
<td><b>Required</b></td>
</tr>
<tr>
<td>[-outseq]<br>(Parameter 2)</td>
<td>Sequence set filename and optional format (output USA)</td>
<td>Writeable sequences</td>
<td><i><*></i>.<i>format</i></td>
</tr>
<tr>
<td>[-dendoutfile]<br>(Parameter 3)</td>
<td>Dendrogram (tree file) from clustalw output file</td>
<td>Output file</td>
<td><i><*></i>.emma</td>
</tr>
<tr bgcolor="#FFFFCC">
<th align="left" colspan=2>Additional (Optional) qualifiers</th>
<th align="left">Allowed values</th>
<th align="left">Default</th>
</tr>
<tr>
<td>-onlydend</td>
<td>Only produce dendrogram file</td>
<td>Toggle value Yes/No</td>
<td>No</td>
</tr>
<tr>
<td>-dend</td>
<td>Do alignment using an old dendrogram</td>
<td>Toggle value Yes/No</td>
<td>No</td>
</tr>
<tr>
<td>-dendfile</td>
<td>Dendrogram (tree file) from clustalw file (optional)</td>
<td>Input file</td>
<td><b>Required</b></td>
</tr>
<tr>
<td>-pwmatrix</td>
<td>The scoring table which describes the similarity of each amino acid to each other.
There are three 'in-built' series of weight matrices offered. Each consists of several matrices which work differently at different evolutionary distances. To see the exact details, read the documentation. Crudely, we store several matrices in memory, spanning the full range of amino acid distance (from almost identical sequences to highly divergent ones). For very similar sequences, it is best to use a strict weight matrix which only gives a high score to identities and the most favoured conservative substitutions. For more divergent sequences, it is appropriate to use 'softer' matrices which give a high score to many other frequent substitutions.
1) BLOSUM (Henikoff). These matrices appear to be the best available for carrying out data base similarity (homology searches). The matrices used are: Blosum80, 62, 45 and 30.
2) PAM (Dayhoff). These have been extremely widely used since the late '70s. We use the PAM 120, 160, 250 and 350 matrices.
3) GONNET . These matrices were derived using almost the same procedure as the Dayhoff one (above) but are much more up to date and are based on a far larger data set. They appear to be more sensitive than the Dayhoff series. We use the GONNET 40, 80, 120, 160, 250 and 350 matrices.
We also supply an identity matrix which gives a score of 1.0 to two identical amino acids and a score of zero otherwise. This matrix is not very useful.</td>
<td><table><tr><td>b</td> <td><i>(blosum)</i></td></tr><tr><td>p</td> <td><i>(pam)</i></td></tr><tr><td>g</td> <td><i>(gonnet)</i></td></tr><tr><td>i</td> <td><i>(id)</i></td></tr><tr><td>o</td> <td><i>(own)</i></td></tr></table></td>
<td>b</td>
</tr>
<tr>
<td>-pwdnamatrix</td>
<td>The scoring table which describes the scores assigned to matches and mismatches (including IUB ambiguity codes).</td>
<td><table><tr><td>i</td> <td><i>(iub)</i></td></tr><tr><td>c</td> <td><i>(clustalw)</i></td></tr><tr><td>o</td> <td><i>(own)</i></td></tr></table></td>
<td>i</td>
</tr>
<tr>
<td>-pairwisedatafile</td>
<td>Comparison matrix file (optional)</td>
<td>Input file</td>
<td><b>Required</b></td>
</tr>
<tr>
<td>-matrix</td>
<td>This gives a menu where you are offered a choice of weight matrices. The default for proteins is the PAM series derived by Gonnet and colleagues. Note, a series is used! The actual matrix that is used depends on how similar the sequences to be aligned at this alignment step are. Different matrices work differently at each evolutionary distance.
There are three 'in-built' series of weight matrices offered. Each consists of several matrices which work differently at different evolutionary distances. To see the exact details, read the documentation. Crudely, we store several matrices in memory, spanning the full range of amino acid distance (from almost identical sequences to highly divergent ones). For very similar sequences, it is best to use a strict weight matrix which only gives a high score to identities and the most favoured conservative substitutions. For more divergent sequences, it is appropriate to use 'softer' matrices which give a high score to many other frequent substitutions.
1) BLOSUM (Henikoff). These matrices appear to be the best available for carrying out data base similarity (homology searches). The matrices used are: Blosum80, 62, 45 and 30.
2) PAM (Dayhoff). These have been extremely widely used since the late '70s. We use the PAM 120, 160, 250 and 350 matrices.
3) GONNET . These matrices were derived using almost the same procedure as the Dayhoff one (above) but are much more up to date and are based on a far larger data set. They appear to be more sensitive than the Dayhoff series. We use the GONNET 40, 80, 120, 160, 250 and 350 matrices.
We also supply an identity matrix which gives a score of 1.0 to two identical amino acids and a score of zero otherwise. This matrix is not very useful. Alternatively, you can read in your own (just one matrix, not a series).</td>
<td><table><tr><td>b</td> <td><i>(blosum)</i></td></tr><tr><td>p</td> <td><i>(pam)</i></td></tr><tr><td>g</td> <td><i>(gonnet)</i></td></tr><tr><td>i</td> <td><i>(id)</i></td></tr><tr><td>o</td> <td><i>(own)</i></td></tr></table></td>
<td>b</td>
</tr>
<tr>
<td>-dnamatrix</td>
<td>This gives a menu where a single matrix (not a series) can be selected.</td>
<td><table><tr><td>i</td> <td><i>(iub)</i></td></tr><tr><td>c</td> <td><i>(clustalw)</i></td></tr><tr><td>o</td> <td><i>(own)</i></td></tr></table></td>
<td>i</td>
</tr>
<tr>