Pairwise sequence alignments Volker Flegel Vassilios Ioannidis VI - - PDF document

VI - 2004 Page 1

Pairwise sequence alignments

Volker Flegel Vassilios Ioannidis

VI - 2004 Page 2

Outline

• Introduction
• Definitions
• Biological context of pairwise alignments
• Computing of pairwise alignments
• Some programs

VI - 2004 Page 3

Importance of pairwise alignments

Sequence analysis tools depending on pairwise comparison

• Multiple alignments
• Profile and HMM making

(used to search for protein families and domains)

• 3D protein structure prediction
• Phylogenetic analysis
• Construction of certain substitution matrices
• Similarity searches in a database

VI - 2004 Page 4

Goal

Sequence comparison through pairwise alignments

• Goal of pairwise comparison is to find conserved regions (if any)

between two sequences

• Extrapolate information about our sequence using the known

characteristics of the other sequence

THIO_EMENI GFVVVDCFATWCGPCKAIAPTVEKFAQTY G ++VD +A WCGPCK IAP +++ A Y ??? GAILVDFWAEWCGPCKMIAPILDEIADEY THIO_EMENI GFVVVDCFATWCGPCKAIAPTVEKFAQTY G ++VD +A WCGPCK IAP +++ A Y ??? GAILVDFWAEWCGPCKMIAPILDEIADEY

THIO_EMENI

SwissProt Extrapolate Extrapolate

???

VI - 2004 Page 5

Do alignments make sense ?

Evolution of sequences

• Sequences evolve through mutation and selection

Selective pressure is different for each residue position in a

protein (i.e. conservation of active site, structure, charge, etc.)

• Modular nature of proteins

Nature keeps re-using domains

• Alignments try to tell the evolutionnary story of the proteins

Relationships

Same Sequence Same 3D Fold Same Origin Same Function VI - 2004 Page 6

Example: An alignment - textual view

• Two similar regions of the Drosophila melanogaster Slit and

Notch proteins

970 980 990 1000 1010 1020 SLIT_DROME FSCQCAPGYTGARCETNIDDCLGEIKCQNNATCIDGVESYKCECQPGFSGEFCDTKIQFC ..:.: :. :.: ...:.: .. : :.. : ::.. . :.: ::..:. :. :. : NOTC_DROME YKCECPRGFYDAHCLSDVDECASN-PCVNEGRCEDGINEFICHCPPGYTGKRCELDIDEC 740 750 760 770 780 790 970 980 990 1000 1010 1020 SLIT_DROME FSCQCAPGYTGARCETNIDDCLGEIKCQNNATCIDGVESYKCECQPGFSGEFCDTKIQFC ..:.: :. :.: ...:.: .. : :.. : ::.. . :.: ::..:. :. :. : NOTC_DROME YKCECPRGFYDAHCLSDVDECASN-PCVNEGRCEDGINEFICHCPPGYTGKRCELDIDEC 740 750 760 770 780 790

VI - 2004 Page 7

Example: An alignment - graphical view

• Comparing the tissue-type and urokinase type plasminogen
• activators. Displayed using a diagonal plot or Dotplot.

Tissue-Type plasminogen Activator Urokinase-Type plasminogen Activator URL: www.isrec.isb-sib.ch/java/dotlet/Dotlet.html

VI - 2004 Page 8

Some definitions

Identity Proportion of pairs of identical residues between two aligned sequences. Generally expressed as a percentage. This value strongly depends on how the two sequences are aligned. Similarity Proportion of pairs of similar residues between two aligned sequences. If two residues are similar is determined by a substitution matrix. This value also depends strongly on how the two sequences are aligned, as well as on the substitution matrix used. Homology Two sequences are homologous if and only if they have a common ancestor. There is no such thing as a level of homology ! (It's either yes or no)

• Homologous sequences do not necessarily serve the same function...
• ... Nor are they always highly similar: structure may be conserved while sequence is not.

VI - 2004 Page 9

More definitions

Consider a set S (say, globins) and a test t that tries to detect members of S (for example, through a pairwise comparison with another globin).

True positive

A protein is a true positive if it belongs to S and is detected by t.

True negative

A protein is a true negative if it does not belong to S and is not detected by t.

False positive

A protein is a false positive if it does not belong to S and is (incorrectly) detected by t.

False negative

A protein is a false negative if it belongs to S and is not detected by t (but should be).

VI - 2004 Page 10 Matches

Definition example

The set of all globins and a test to identify them

True positives True negatives False positives False negatives

Consider:

• a set S (say, globins: G)
• a test t that tries to detect members of S

(for example, through a pairwise comparison with another globin). Globins

G G G G G G G G X X X X X

VI - 2004 Page 11

Even more definitions

Sensitivity

Ability of a method to detect positives, irrespective of how many false positives are reported.

Selectivity

Ability of a method to reject negatives, irrespective of how many false negatives are rejected.

True positives True negatives False positives False negatives Greater sensitivity Less selectivity Less sensitivity Greater selectivity VI - 2004 Page 12

Pairwise sequence alignment

Concept of a sequence alignment

• Pairwise Alignment:

Explicit mapping between the residues of 2 sequences

– Tolerant to errors (mismatches, insertion / deletions or indels) – Evaluation of the alignment in a biological concept (significance)

Seq AGARFIELDTHELASTFA-TCAT ||||||||||| || |||| Seq BGARFIELDTHEVERYFASTCAT Seq AGARFIELDTHELASTFA-TCAT ||||||||||| || |||| Seq BGARFIELDTHEVERYFASTCAT

errors / mismatches insertion deletion

VI - 2004 Page 13

Pairwise sequence alignement

Number of alignments

• There are many ways to align two sequences
• Consider the sequence fragments below: a simple alignment

shows some conserved portions but also:

CGATGCAGACGTCA |||||||| CGATGCAAGACGTCA CGATGCAGACGTCA |||||||| CGATGCAAGACGTCA CGATGCAGACGTCA |||||||| CGATGCAAGACGTCA CGATGCAGACGTCA |||||||| CGATGCAAGACGTCA

• Number of possible alignments for 2 sequences of length 1000 residues:

more than 10600 gapped alignments

(Avogadro 1024, estimated number of atoms in the universe 1080) VI - 2004 Page 14

Alignement evaluation

What is a good alignment ?

• We need a way to evaluate the biological meaning of a given alignment
• Intuitively we "know" that the following alignment:

is better than:

CGAGGCACAACGTCA ||| ||| |||||| CGATGCAAGACGTCA CGAGGCACAACGTCA ||| ||| |||||| CGATGCAAGACGTCA ATTGGACAGCAATCAGG | || | | ACGATGCAAGACGTCAG ATTGGACAGCAATCAGG | || | | ACGATGCAAGACGTCAG

• We can express this notion more rigorously, by using a

scoring system

VI - 2004 Page 15

Scoring system

Simple alignment scores

• A simple way (but not the best) to score an alignment is to count 1 for

each match and 0 for each mismatch.

Score: 12

CGAGGCACAACGTCA ||| ||| |||||| CGATGCAAGACGTCA CGAGGCACAACGTCA ||| ||| |||||| CGATGCAAGACGTCA ATTGGACAGCAATCAGG | || | | ACGATGCAAGACGTCAG ATTGGACAGCAATCAGG | || | | ACGATGCAAGACGTCAG

Score: 5

VI - 2004 Page 16

Introducing biological information

Importance of the scoring system discrimination of significant biological alignments

• Based on physico-chemical properties of amino-acids

Hydrophobicity, acid / base, sterical properties, ... Scoring system scales are arbitrary

• Based on biological sequence information

Substitutions observed in structural or evolutionary alignments of

well studied protein families

Scoring systems have a probabilistic foundation Substitution matrices

• In proteins some mismatches are more acceptable than others
• Substitution matrices give a score for each substitution of one amino-

acid by another

VI - 2004 Page 17

Substitution matrices (log-odds matrices)

Example matrix

PAM250

From:

• A. D. Baxevanis, "Bioinformatics"

(Leu, Ile): 2 (Leu, Cys): -6 ...

• Positive score: the amino acids are

similar, mutations from one into the other occur

more often then expected by chance during evolution

• Negative score: the amino acids are

dissimilar, the mutation from one into the other

• ccurs less often then expected by chance during

evolution

• chance

by expected

• bserved

log

• For a set of well known proteins:
• Align the sequences
• Count the mutations at each position
• For each substitution set the score to the

log-odd ratio VI - 2004 Page 18

Matrix choice

Different kind of matrices

• PAM series

(Dayhoff M., 1968, 1972, 1978)

Percent Accepted Mutation. A unit introduced by Dayhoff et al. to quantify the amount of evolutionary change in a protein sequence. 1.0 PAM unit, is the amount of evolution which will change, on average, 1% of amino acids in a protein sequence. A PAM(x) substitution matrix is a look-up table in which scores for each amino acid substitution have been calculated based on the frequency of that substitution in closely related proteins that have experienced a certain amount (x) of evolutionary divergence.

Based on 1572 protein sequences from 71 families Old standard matrix:

PAM250

VI - 2004 Page 19

Matrix choice

Different kind of matrices

• BLOSUM series

(Henikoff S. & Henikoff JG., PNAS, 1992)

Blocks Substitution Matrix. A substitution matrix in which scores for each position are derived from

• bservations of the frequencies of substitutions in blocks of local alignments in

related proteins. Each matrix is tailored to a particular evolutionary distance. In the BLOSUM62 matrix, for example, the alignment from which scores were derived was created using sequences sharing no more than 62% identity. Sequences more identical than 62% are represented by a single sequence in the alignment so as to avoid over-weighting closely related family members.

Based on alignments in the BLOCKS database Standard matrix:

BLOSUM62

VI - 2004 Page 20

Matrix choice

Limitations

• Substitution matrices do not take into account long range

interactions between residues.

• They assume that identical residues are equal (whereas in reallife

a residue at the active site has other evolutionary constraints than the same residue outside of the active site)

• They assume evolution rate to be constant.

VI - 2004 Page 21

Alignment score

Amino acid substitution matrices

• Example:

PAM250

• Most used:

Blosum62

Raw score of an alignment

TPEA _| | APGA TPEA _| | APGA

Score = 1 = 9 + 6 + + 2

VI - 2004 Page 22

Gaps

Insertions or deletions

• Proteins often contain regions where residues have been inserted or

deleted during evolution

• There are constraints on where these insertions and deletions can

happen (between structural or functional elements like: alpha helices, active site, etc.)

Gaps in alignments

GCATGCATGCAACTGCAT ||||||||| GCATGCATGGGCAACTGCAT GCATGCATGCAACTGCAT ||||||||| GCATGCATGGGCAACTGCAT

can be improved by inserting a gap

GCATGCATG--CAACTGCAT ||||||||| ||||||||| GCATGCATGGGCAACTGCAT GCATGCATG--CAACTGCAT ||||||||| ||||||||| GCATGCATGGGCAACTGCAT

VI - 2004 Page 23

Gap opening and extension penalties

Costs of gaps in alignments

• We want to simulate as closely as possible the evolutionary

mechanisms involved in gap occurence.

Example

• Two alignments with identical number of gaps but very different gap
• distribution. We may prefer one large gap to several small ones

(e.g. poorly conserved loops between well-conserved helices)

CGATGCAGCAGCAGCATCG |||||| ||||||| CGATGC------AGCATCG CGATGCAGCAGCAGCATCG |||||| ||||||| CGATGC------AGCATCG CGATGCAGCAGCAGCATCG || || |||| || || | CG-TG-AGCA-CA--AT-G CGATGCAGCAGCAGCATCG || || |||| || || | CG-TG-AGCA-CA--AT-G

gap opening

Gap opening penalty

• Counted each time a gap is opened in an alignment

(some programs include the first extension into this penalty) gap extension

Gap extension penalty

• Counted for each extension of a gap in an alignment

VI - 2004 Page 24

Gap opening and extension penalties

Example

• With a match score of 1 and a mismatch score of 0
• With an opening penalty of 10 and extension penalty of 1, we have

the following score:

CGATGCAGCAGCAGCATCG |||||| ||||||| CGATGC------AGCATCG CGATGCAGCAGCAGCATCG |||||| ||||||| CGATGC------AGCATCG CGATGCAGCAGCAGCATCG || || |||| || || | CG-TG-AGCA-CA--AT-G CGATGCAGCAGCAGCATCG || || |||| || || | CG-TG-AGCA-CA--AT-G

gap opening

13 x 1 - 10 - 6 x 1 = -3

gap extension

13 x 1 - 5 x 10 - 6 x 1 = -43

VI - 2004 Page 25

Statistical evaluation of results

Alignments are evaluated according to their score

• Raw score

It's the sum of the amino acid substitution scores and gap

penalties (gap opening and gap extension)

Depends on the scoring system (substitution matrix, etc.) Different alignments should not be compared based only on

the raw score

• It is possible that a "bad" long alignment gets a better raw score than a very good

short alignment.

We need a normalised score to compare alignments ! We need to evaluate the biological meaning of the score (p-value, e-value).

• Normalised score

Is independent of the scoring system Allows the comparison of different alignments Units: expressed in bits

VI - 2004 Page 26

Statistical evaluation of results

Distribution of alignment scores - Extreme Value Distribution

• Random sequences and alignment scores

Sequence alignment scores between random sequences are

distributed following an extreme value distribution (EVD).

Ala Val ... Tr p

score

• bs

score x score y ... ...

Ala Val ... Tr p Random sequences Pairwise alignments Score distribution

VI - 2004 Page 27

score y: our alignment is very improbable to

• btain with random

sequences

Statistical evaluation of results

Distribution of alignment scores - Extreme Value Distribution

• High scoring random alignments have a low probability.
• The EVD allows us to compute the probability with which our

biological alignment could be due to randomness (to chance).

• Caveat: finding the threshold of significant alignments.

score

score x: our alignment has a great probability

• f being the result of

random sequence similarity

Threshold significant alignment

VI - 2004 Page 28

Statistical evaluation of results

Statistics derived from the scores

• p-value

Probability that an alignment with this score occurs by chance

in a database of this size

The closer the p-value is towards 0, the better the alignment

• e-value

Number of matches with this score one can expect to find by

chance in a database of this size

The closer the e-value is towards 0, the better the alignment

• Relationship between e-value and p-value:

In a database containing N sequences

e = p x N

100% 0% N

VI - 2004 Page 29

Diagonal plots or Dotplot

Concept of a Dotplot

• Produces a graphical representation of similarity regions.
• The horizontal and vertical dimensions correspond to the

compared sequences.

• A region of similarity stands out as a diagonal.

Tissue-Type plasminogen Activator Urokinase-Type plasminogen Activator

VI - 2004 Page 30

Dotplot construction

Simple example

• A dot is placed at each position where two residues match.

The colour of the dot can be chosen according to the substitution

value in the substitution matrix

THEFA-TCAT ||||| |||| THEFASTCAT THEFA-TCAT ||||| |||| THEFASTCAT

Note

• This method produces dotplots with too much noise to be useful

The noise can be reduced by calculating a score using a window

• f residues

The score is compared to a threshold or stringency

VI - 2004 Page 31

Dotplot construction

Window example

• Each window of the first sequence is aligned (without gaps) to each

window of the 2nd sequence

• A colour is set into a rectangular array according to the score of the

aligned windows

THE ||| THE THE ||| THE

Score: 23

THE HEF THE HEF

Score: -5

CAT THE CAT THE

Score: -4

HEF THE HEF THE

Score: -5 VI - 2004 Page 32

Dotplot limitations

It's a visual aid.

The human eye can rapidly identify similar regions in sequences.

It's a good way to explore sequence organisation. It does not provide an alignment.

Tissue-Type plasminogen Activator Urokinase-Type plasminogen Activator

VI - 2004 Page 33

Relationship between alignment and dotplot

• An alignment can be seen as a path through the dotplot diagramm.

Creating an alignment

Seq B A-CA-CA | || | Seq A ACCAAC- Seq B A-CA-CA | || | Seq A ACCAAC- Seq B ACA--CA | Seq A A-CCAAC Seq B ACA--CA | Seq A A-CCAAC VI - 2004 Page 34

Finding an alignment

Alignment algorithms

• An alignment program tries to find the best alignment between two

sequences given the scoring system.

• This can be seen as trying to find a path through the dotplot diagram including all (or the

most visible) diagonals.

Alignement types

• Global

Alignment between the complete sequence A and the complete sequence B

• Local

Alignment between a sub-sequence of A an a sub- sequence of B

Computer implementation (Algorithms)

• Dynamic programing
• Global

Needleman-Wunsch

• Local

Smith-Waterman

VI - 2004 Page 35

Global alignment (Needleman-Wunsch)

Example Global alignments are very sensitive to gap penalties Global alignments do not take into account the modular nature of

proteins

Tissue-Type plasminogen Activator Urokinase-Type plasminogen Activator

Global alignment: VI - 2004 Page 36

Local alignment (Smith-Waterman)

Example Local alignments are more sensitive to the modular nature of

proteins

They can be used to search databases

Tissue-Type plasminogen Activator Urokinase-Type plasminogen Activator

Local alignments:

VI - 2004 Page 37

Optimal alignment extension

How to extend optimaly an optimal alignment

• An optimal alignment up to positions i and j can be extended in 3 ways.
• Keeping the best of the 3 guarantees an extended optimal alignment.

Seq A a1 a2 a3 ... ai-1 ai Seq B b1 b2 b3 ... bj-1 bj Seq A a1 a2 a3 ... ai-1 ai Seq B b1 b2 b3 ... bj-1 bj

• We have the optimal alignment extended from i and j by one residue.

Seq A a1 a2 a3 ... ai-1 ai Seq B b1 b2 b3 ... bj-1 bj Seq A a1 a2 a3 ... ai-1 ai Seq B b1 b2 b3 ... bj-1 bj ai+1 bj+1 ai+1 bj+1

Score = Scoreij + Substi+1j+1

Seq A a1 a2 a3 ... ai-1 ai Seq B b1 b2 b3 ... bj-1 bj Seq A a1 a2 a3 ... ai-1 ai Seq B b1 b2 b3 ... bj-1 bj ai+1

• ai+1
• Score = Scoreij - gap

Seq A a1 a2 a3 ... ai-1 ai Seq B b1 b2 b3 ... bj-1 bj Seq A a1 a2 a3 ... ai-1 ai Seq B b1 b2 b3 ... bj-1 bj

• bj+1
• bj+1

Score = Scoreij - gap VI - 2004 Page 38

Exact algorithms

Simple example (Needleman-Wunsch)

• Scoring system:

Match score:

2

Mismatch score:

• 1

Gap penalty:

• 2

Note

• We have to keep track of the origin of the score for each element in

the matrix. This allows to build the alignment by traceback when the matrix has

been completely filled out.

• Computation time is proportional to the size of sequences (n x m).

20

• 4
• 4

0 - 2 0 - 2 2 + 2

20

• 4
• 4

F(i-

1,j)

F(i,j)

s(xi,yj)

F(i-1,j-

1)

• d

F(i,j-

1)

• d

F(i,j): score at position i, j s(xi,yj): match or mismatch score (or substitution matrix value) for residues xi and yj d: gap penalty (positive value)

GA-TTA || || GAATTC GA-TTA || || GAATTC

VI - 2004 Page 39

Algorithms for pairwise alignments

Web resources

• LALIGN - pairwise sequence alignment:

www.ch.embnet.org/software/LALIGN_form.html

• PRSS - alignment score evaluation:

www.ch.embnet.org/software/PRSS_form.html

Concluding remarks

• Substitution matrices and gap penalties introduce biological

information into the alignment algorithms.

• It is not because two sequences can be aligned that they share

a common biological history. The relevance of the alignment must be assessed with a statistical score.

• There are many ways to align two sequences.

Do not blindly trust your alignment to be the only truth. Especially gapped regions may be quite variable.

• Sequences sharing less than 20% similarity are difficult to align:

You enter the Twilight Zone (Doolittle, 1986) Alignments may appear plausible to the eye but are no longer

statistically significant.

Other methods are needed to explore these sequences (i.e:

profiles)