Sequence Alignment Methods And Algorithms

Algorithms for Sequence Alignments A. Brüngger, Labhead Bioinformatics Novartis Pharma AG [email protected]

Algorithms for Sequence Alignments • Definition of Sequence Alignment and Origins of Sequence Similarity – – – –

Global Alignment Local Alignment Alignment of Pairs of Sequences Alignment of Multiple Sequences

• Methods for Pairwise Sequence Alignments – Dot matrix – Scoring Matrices – Dynamic Programming: Smith-Waterman, Needleman-Wunsch

• Methods for Multiple Sequence Alignment – Dynamic Programming – Progressive multiple alignments: CLUSTALW, PILEUP

• Database Searching for Similar Sequences – FASTA, BLAST – Patscan – Hidden Markov Models

Outline follows: David W. Mount, "Bioionformatics - Sequence and Genome Analysis“ Cold Spring Harbour Laboratory Press, 2001. Online: http://www.bioinformaticsonline.or g

Rational for Sequence Analysis, Origins of Sequence Similarity Function Structure

Similar sequence leads to similar function

Protein

Sequence Analysis as the basic tool to discover functional, structural,evolutionary information in biological sequences

DNA

Sequence A

y Steps

Sequence B

x Steps

common ancestor sequence

Evolutionary relationship between two similar sequences and a possible common ancestor. The number of steps to convert one sequence into the other is the "evolutionary" distance between the sequences (x + y). Usually, the ancestor sequence is not available, only (x + y) can be computed.

Origins of Homology  Significance of Sequence Alignments

Possible Origins of Sequence Homology: • • • •

orthologs (panel A and B) a1 in species I and a1 in species II (same ancestor!) paralogs (panel A and B) a1 and a2 (arose from gene duplication event) analogs (panel C): different genes converge to same function by different evolutionary paths transfer of genetic material (panel D) between different species

Homology vs. Similarity • •

Similarity can be computed (by sequence alignments) Homology is deduced (e.g. from similarity, but also from other evidence!)

Definition of Sequence Alignment Computational procedure (“algorithm”) for comparing two/many sequences • identify series of identical residues or patterns of identical residues that appear in the same order in the sequences • visualized by writing sequences as follows:

MLGPSSKQTGKGS-SRIWDN* || | ||| | | MLN-ITKSAGKGAIMRLGDA*

Pairwise Global Alignment (over whole length of sequences)

GKG ||| GKG

Pairwise Local Alignment (similar parts of sequences)

• sequence alignment is an optimiztion problem bringing as many identical residues as possible into corresponding positions

Pairwise Sequence Alignment Alignment Score: Quantify the Quality of an alignment • simple scoring system (eg. for DNA sequences) – when two identical residues are aligned: +1 – when two non-identical residues are aligned: -1 – when a gap is introduced -1

• total score of alignment is the sum of the scores in each position:

agaag-tagattcta || || ||| || || aggaggtag-tt-ta

Compute Score:

•11 matches •1 mismatch •3 gaps Score = 11 - 1 -3 = 7

Alignment is 'optimal' when assigned score is maximal Note: Scoring scheme defines 'optimal' alignment

Algorithms for Pairwise Sequence Alignments: Dot Plot Algorithm introduced by Gibbs and McIntyre (1970): • write one sequence on top from left to right • write other sequence on the left from top to bottom • place a dot whenever two residues are identical; variants use colors • visual inspection (variants: identify diagonals with “many” dots

a g a c c t a g a * * * g * * a * * * c * * c * * t *

Algorithms for Pairwise Sequence Alignments: Dot Plot Example Dot Plots of DNA and Protein of Phage l cI (horizontal) and P22 c3 (vertical)

DNA

Protein

Removing “noise”: Plot a dot only if 7 ("stringency") out of the next 11 ("window size") residues are identical

Algorithms for Pairwise Sequence Alignments: Dot Plot Example: Human LDL receptor against itself.

Window size 1, stringency 1.

Window size 23, stringency 7.

"Repeats" in first part of sequence.

Algorithms for Pairwise Sequence Alignments: Dot Plot • Strengths: – – – –

visualization of sequence similarity finding direct and inverted repeats in sequences finding self-complementary regions in RNA (secondary structures) simple to compute, simple to visualize

• Implementations DNA Strider (Macintosh) DOTTER (UNIX X-Windows) GCG "COMPARE", "DOTPLOT" online SIB: http://www.isrec.isb-sib.ch/java/dotlet/Dotlet.html

• Computational Complexity: – two sequences of length n, m: O(n·m)

compute time

– – – –

sequence length

Scoring Matrices for Proteins Observation: Protein (Function, Structure) is preserved when substitutions in certain AAs happen • scoring system less simple in case of proteins • Example: The same protein in three different species

A W T V A S A V R T S I A Y T V A S A V R T S I A W T V A A A V L T S I A

B W to Y L to R A to S

Organism A Organism B Organism C

C Therefore: Assigning L to R scores higher than assigning it to another AA

Scoring Matrices for Proteins Generalization: Dayhoff PAM Weight Matrices (percent accepted mutation) – observed mutations during a unit of evolutionary time (1 PAM ~ time that is required that an AA is changed with probability 1%) – initial PAM matrix ("PAM1"): derived from 1572 changes in 71 groups of protein sequences that are at least 85% similar – observed in proteins that have the same function => "accepted" mutations

C S . .

C S T P .. fcc fcs fct fcp fsc fss fst fsp

Scoring Matrices for Proteins Extrapolation to longer evolutionary distances possible!

M

R

K

H

PAM1: p(HM) = 0 PAM2: p(HM) = p(HK)*p(KM) +p(HR)*p(RM)

Scoring Matrices for Proteins • For each pair of residues, a score is given by a "scoring matrix" • Scoring of two pairs of residues depends on: – frequency of the two residues – exchange that tends to preserves function should have high score

• Scoring matrices are constructed by empirically evaluating (manually constructed) alignments of closely related proteins • PAM – introduced by Margarete Dayhoff, 1978 – inspecting 1572 changes in 71 groups of protein sequences – different matrices for varying degree of similarity (or, evolutionary distance)

• BLOSUM – inspection of about 2'000 conserved AA-stretches, "BLOCKS" Sequence 1: Sequence 2: score

V V 4

D E 2

S S L 4 -2

C C 4

Y Y 4

Total 16

Dynamic Programming Approach to Sequence Alignments • Description by Example: – input:

two AA sequences “VDSCY” and “VESLCY” scoring “matrix”: match +4, mismatch +2, gap -2 – Some possible alignments and their score: VDSCY VD-SCY VDS-CY | | | | | | || VESLCY V-ESLCY VESLCY 14 8 16 – Observation: Longer alignments can be derived from shorter new score VDS-CY VESLCY 16 VDS-C VESLC 12

= old score

=

VDS-C VESLC 12

=

VDSVESL 8

+ score of new pair

+

Y Y 4

+

C C 4

Dynamic Programming Approach to Sequence Alignments • Alignment: Path in matrix from "top left" to "bottom right" V D S C Y Seq 1: Seq 2:

V D S - C Y V E S L C Y

V E S L C Y possible extensions of "current" alignment ( )

• Three situations for extending previous alignment: introduce gap in sequence 1

introduce gap in sequence 2

Dynamic Programming Approach to Sequence Alignments

• In each matrix element, we can write the score of the best alignment up to this element, and a pointer to the shorter alignment it is derive from V D S C Y

Seq 1: Seq 2:

V V 4

D E 2

S S L 4 -2

C C 4

Y Y 4

match +4, mismatch +2, gap -2

V D S C Y V E S L C Y

4 2 2 2 2 2 2 2 2 2

V E S L C Y

4

Goal is to compute this path!

6 10 12 16

V D S C Y

two possibilities: 4 + (E->S) + gap = 4 2 + (E->S) = 4 choose one (first)

three possibilities: 4 + (L->D) + 2 gaps = 2 2 + (L->D) + 1 gap = 2 2 + (L->D) = 4 choose third (score max)

fill first row/column

V E S L C Y

4 2 2 2 2 2 6 4 4 4 2 4 2 4 2 4 2 4

Dynamic Programming Approach to Sequence Alignments • Formally:

Si-1, j-1 + s(ai -> bi) Sij = max

max [x≥1] (Si-x,j - wx) max [y≥1] (Si,j-y - wy)

– – – –

ai, bj : Residue of Sequence A at position i and Residue of Sequence B at position j wx : Penalty for introducing gap of length x Sij : Score of alignment at position (i, j) s(a->b) : Score for aligning a and b (known from scoring matrix)

• Iterate process until whole matrix is filled with scores and backpointers • Choose maximum score in last column or row • Follow pointers to construct alignment

Dynamic Programming Approach to Sequence Alignments • Global Alignment Needleman and Wunsch (1970) • Local Alignment: Smith-Waterman (minor modification)

Dynamic Programming Approach to Sequence Alignments • Implementations available on the web

Dynamic Programming Approach to Sequence Alignments • Strengths: – – – – –

finds optimal (mathematically best) alignment suited for both, local and global alignments global alignment: Needleman-Wunsch local alignment: Smith-Waterman statistical significance can be attached (recompute alignment when one or both sequences are randomly changed)

• Implementations – LALIGN – GCG "GAP" (global) and BESTFIT (local) – ...

• Complexity: – two sequences of length n, m – time: O(n·m), space: O(n·m) – space is crucial • example: matrix element 4 bytes, n=m=10000, space requirement: 400MB • can be improved: trade time for space

Multiple Sequence Alignment • When aligning k (k>2) sequences, the dynamic programming idea can theoretically be used: – instead of a two-dimensional scoring schema, a k-dimensional one is used – instead of two-dimensional residue substitution matrix, a kdimensional one is used

• However, this is not feasible in practical terms – size of these matrices increase too rapidly! – example: 5 sequences, 100 residues each, matrix size = 100^5 = 10^10 in terms of computer memory: 10 GB

• Consequence: Optimal alignment for multiple sequences is not computable! • Approximative methods used (“heuristics”)

Multiple Sequence Alignment: Progressive methodes

• Idea for progressive methodes: – – – –

step step step step

1 2 3 4

: : : :

Align two of the sequences Compute "consensus" sequence Align a third sequence to the consensus Repeat steps 2 and 3

Multiple Sequence Alignment: Progressive methodes, CLUSTALW • Basic steps – – – – –

step 1 : Perform pairwise sequence alignments of all possible pairs step 2 : Use scores to produce a “phylogenetic tree” step 3 : align the sequences sequentially, guided by the tree the most closely related sequences are aligned first other sequences or groups of sequences are added

• Example: Phylogenetic tree for 4 sequences A B C D A 100 90 85 90 B 100 95 80 C 100 97 D 100 Percentage Identities in pairwise alignment

A

90

B

90 85

80 95

C

97

D

Join sequences with highest similarity

Multiple Sequence Alignment: Progressive methodes, CLUSTALW • Computation of Phylogenetic tree Join sequences with highest similarity

A B C D A 100 95 80 90 B 100 95 85 C 100 97 D 100 Percentage Identities in pairwise alignment

A

95

A

B

95

B

90 80

85

(80+90)/2 = 85

95 C

97

(85+95)/2 = 90

D C, D

Join sequences with highest similarity

A, B

A B C D

(85+90)/2 = 87.5

C, D

Multiple Sequence Alignment: Progressive methodes, CLUSTALW • Example: Seven Globins from SWISSPROT

Multiple Sequence Alignment: Progressive methodes, CLUSTALW • Available Implementations

Riddle: Probability that a short string occurs in a longer text

• Given:

an alphabet A = {a, c, g, t} a random string s of length k over A a random text t of length n (n>k) over A

• Find probability p, by which the string s occurs in the text t. • Find expected number of occurrences e of s in t.

• Example: s = tggtacaaatgttct (glucocorticoid response element GR t = 10’000 bp (‘promoter’, upstream DNA to start codon

Basics of Sequence DB Searches: Definition

• Given:

one “short” sequence q (query) sequence DB, conceptually one long sequence s

(subject) • Find occurrences of “similar” sequences to q in s. • Attach to each similar occurrence a statistical significance. • Example: s = human genomic DNA, 3·109 b q = tggtacaaatgttct (glucocorticoid response element GRE) • Dynamic programming approach not feasible!

Basics of Sequence DB Searches: Detection of identical kmers • Idea:

identify identical k-mers in s and q (“seed”) expand alignment from seed in both directions

• Example: q

MAAARLCLSLLLLSTCVALLLQPLLGAQGAPLEPVYPGDNATPEQMAQYAADLRRYINMLTRPRYGKRHKEDTLAFSEWGS || | |||| | || |||||||| | |||||| |||| ||||||||| |||||| ||||||||| | ... MAVAYCCLSLFLVSTWVALLLQPLQGTWGAPLEPMYPGDYATPEQMAQYETQLRRYINTLTRPRYGKRAEEENTGGLP...

1 2 Gonnet120 score: 412, 76 % identities

3 no more increase in score

• BLAST • HSP, “high scoring pair” • gapped alignment • starting extension also from similar (and not only identical) seeds

Basics of Sequence DB Dearches: Detection of identical kmers • Precompute position of all k-mers in DB sequence • Indexing all Peptides of length k in "database“ 0 1 2 3 • Example: 1234567890123456789012345678901234 MAAARLCLSLLLLSTCVALLLQPLLGAQGAPLEP MAAAR AAARL AARLC .... APLEP Sorted: AAARL 2 AARLC 3 APLEP 30 ... MAAAR 1 ... VALLL 17

• For each peptide of length k in "query", the identical peptides in "database" are identified efficiently (binary search in sorted list) • For identical pairs, extension step in both directions is performed

Basics of Sequence DB Dearches: Detection of identical kmers • Indexing all Peptides of length k in "database“: some refinements 0 1 2 3 1234567890123456789012345678901234 MAAARLCVALLLLSTCVALLLQPLLGAQGAPLEP 8

Pointers to previous occurrences

List of all words of length k and last occurrence in query: Simple, yet efficient data-structure: AAAAA - array of integers (size=length of db) AAAAC AAAAD - array of integers (size=number of words with length k .... AARLC 3 Book-keeping: .... - more than one db/query sequence APLEP 30 - build database chunks that fit into main memory ... (speeds up computation 1000x) ... VALLL 17 ... Extension step: optimizations

•

For each peptide of length k in "query", the position in the wordlist can be easily computed (no binary search!)

Significance of matches: DNA case

• issue: searching with short query vs. large database  found mat could have occurred by pure chance • assume equal distribution of c,g,a,t • what is ...

– the probability q, that sequence B (len=m) is contained in sequence A (len=n)? – the expected length of a common subsequence of two sequences? – the expected score when locally aligning two sequences of length n, m – Solution for first: • a. There are (n-m) 'words' of length m in sequence A • b. In total, there are 4^m sequences of length m • c. p = (n-m) / 4^m

This is wrong. Why?

Not independent!

Significance of matches: Statistics

• statistics – the statistical distribution of alignment scores found in a DB search follows the extreme value distribution (not normal distribution) – extreme value distribution changes with length of sequences and their residual composition – scores of actual database search results are plotted vs. expected scores (FASTA) – BLAST computes E-Value (number of expected hits with this score, when comparing the query with unrelated database sequences)

Putting it together: BLAST2

A W T V A S A V R T S I

• (optional) filtering for low complexity region in query • all words of length 3 are listed: AWT VAS AVR TSI | WTV ASA VRT | TVA SAV RTS

• to each word, ~50 'high scoring' additional words are added: AWT AWA TWA ART ...

VAS IAS LAS ITS ...

AVR TVR AIR AVS ...

TSI | WTV ASA VRT | TVA SAV RTS ... ... ...

• matching words are found in DB (with datastructure as described before) • ungapped alignment constructed from word matches, 'HSP' • statistics determines, whether HSP is significant • SW-alignment for significant HSPs

An example: comparing mus chr X vs. hum chr X, syntenic regions mouse chromosome X (147Mbp)

human chromosome X (149Mbp)

Each dot: conserved stretch of AA HSP, high scoring pair

Conclusions

• • • • • • •

types of sequence alignments: pairwise, multiple, query vs. database local and global alignment scoring matrices: attach score to each pair of residues (allow „similar“ residues to be aligned) sequence alignment problem is mathematical optimization problem: find optimal alignment: maximise score optimal alignment in practical terms only feasible for pairwise alignment Heuristics for multiple sequence alignments: progressive alignment guided by all pairwise alignments Sequence database searches: – efficiently find occurrence of query k-mers in db sequences (seeds) – expand (ungapped) HSP from seeds – attach statistical significance