Transcript Slides - Fenyo Lab

Previous Lecture: Probability
This Lecture
Introduction to Biostatistics and Bioinformatics
Sequence Alignment Concepts
Sequence Alignment
Stuart M. Brown, Ph.D.
Center for Health Informatics and Bioinformatics
NYU School of Medicine
Slides/images/text/examples borrowed liberally from:
Torgeir R. Hvidsten,
Michael Schatz,
Bill Pearson,
Fourie Joubert, others …
Learning Objectives
• Identity, similarity, homology
• Analyze sequence similarity by dotplots
• window/stringency
Alignment of text strings by edit distance
Scoring of aligned amino acids
Gap penalties
Global vs. local alignment
Dynamic Programming (Smith
Waterman)
• FASTA method
•
•
•
•
•
Why Compare Sequences?
• Identify sequences found in lab experiments
• What is this thing I just found?
• Compare new genes to known ones
• Compare genes from different species
• information about evolution
• Guess functions for entire genomes full of
new gene sequences
• Map sequence reads to a Reference Genome
(ChIP-seq, RNA-seq, etc.)
Are there other sequences like
this one?
1) Huge public databases - GenBank, Swissprot, etc.
2) “Sequence comparison is the most powerful and
reliable method to determine evolutionary
relationships between genes” -R. Pearson
3) Similarity searching is based on alignment
4) BLAST and FASTA provide rapid similarity
searching
a. rapid = approximate (heuristic)
b. false + and - scores
Similarity ≠ Homology
1) 25% similarity ≥ 100 AAs is strong
evidence for homology
2) Homology is an evolutionary statement
which means “descent from a common
ancestor”
– common 3D structure
– usually common function
– homology is all or nothing, you cannot say
"50% homologous"
Similarity is Based on Dot Plots
1) two sequences on vertical and horizontal
axes of graph
2) put dots wherever there is a match
3) diagonal line is region of identity
(local alignment)
4) apply a window filter - look at a group of
bases, must meet % identity to get a dot
Simple Dot Plot
G A T CA A C TG A C GT A
G
T
T
C
A
G
C
T
G
C
G
T
A
C
Window / Stringency
Score = 11
PTHPLASKTQILPEDLASEDLTI
PTHPLAGERAIGLARLAEEDFGM

Filtering
Score = 11
PTHPLASKTQILPEDLASEDLTI
PTHPLAGERAIGLARLAEEDFGM
Score = 7
PTHPLASKTQILPEDLASEDLTI
PTHPLAGERAIGLARLAEEDFGM

Window = 12
Stringency = 9
Dot plot filtered with 4 base
window and 75% identity
G A T CA A C TG A C GT A
G
T
T
C
A
G
C
T
G
C
G
T
A
C
Dot plot of real data
Dotplot
(Window = 130 / Stringency = 9)
Hemoglobin
-chain
Hemoglobin -chain
Dotplot
(Window = 18 / Stringency = 10)
Hemoglobin
-chain
Hemoglobin -chain
How to Align Sequences?
• Manually line them up and count?
• an alignment program can do it for you
• or a just use a text editor
GATGCCATAGAGCTGTAGTCGTACCCT <—
—> CTAGAGAGC-GTAGTCAGAGTGTCTTTGAGTTCC
• Dot Plot
– shows regions of similarity as diagonals
Percent Sequence Identity
• The extent to which two nucleotide or amino
acid sequences are invariant
AC C TG A G – AG
AC G TG – G C AG
mismatch
indel
70% identical
Hamming Distance
• The minimum number of base changes that will
convert one (ungapped) sequence into another
Python function hamming_distance
def hamming_distance(s1, s2):
#Return the Hamming distance between equal-length sequences
if len(s1) != len(s2):
raise ValueError("Undefined for sequences of unequal length")
return sum(ch1 != ch2 for ch1, ch2 in zip(s1, s2))
• The Hamming distance is named after Richard
Hamming, who introduced it in his fundamental
paper on Hamming codes: “Error detecting and error
correcting codes” (1950) Bell System Technical Journal 29 (2): 147–
160.
Hamming Dist can be unrealistic
v: ATATATAT
w: TATATATA
Hamming Dist = 8 (no gaps, no shifts)
• Levenshtein (1966) introduced edit distance
v = _ATATATAT
w = TATATATA_
edit distance: d(v, w) = 2
Levenshtein, Vladimir I. (February 1966). "Binary codes capable of correcting
deletions, insertions, and reversals". Soviet Physics Doklady 10 (8): 707–710.
Affine Gap Penalties
Gap Penalites
• With unlimited gaps (no penalty), unrelated
sequences can align (especially DNA)
• Gap should cost much more than a mismatch
• Multi-base gap should cost only a little bit more
than a single base gap
• Adding an additional gap near another gap
should cost more (not implemented in most algorithms)
• Score for a gap of length x is: -(p + σx)
• p is gap open penalty
• σ is gap extend penalty
Global vs Local similarity
1) Global similarity uses complete aligned sequences
- total % matches
- Needleman & Wunch algorithm
2) Local similarity looks for best internal matching
region between 2 sequences
- find a diagonal region on the dotplot
– Smith-Waterman algorithm
– BLAST and FASTA
3) dynamic programming
– optimal computer solution, not approximate
Global vs. Local Alignments
Smith-Waterman Method
[Essentially finding the diagonals on the dotplot]
•Basic principles of dynamic programming
- Creation of an alignment path matrix
- Stepwise calculation of score values
- Backtracking (evaluation of the optimal path)
Creation of an alignment path matrix
Idea: Build up an optimal alignment using previous
solutions for optimal alignments of smaller
subsequences
• Construct matrix F indexed by i and j (one index for each
sequence)
• F(i,j) is the score of the best alignment between the initial
segment x1...i of x up to xi and the initial segment y1...j of y up
to yj
• Build F(i,j) recursively beginning with F(0,0) = 0
Michael Schatz
Creation of an alignment path matrix
• If F(i-1,j-1), F(i-1,j) and F(i,j-1) are known we can calculate
F(i,j)
• Three possibilities:
• xi and yj are aligned, F(i,j) = F(i-1,j-1) + s(xi ,yj)
• xi is aligned to a gap, F(i,j) = F(i-1,j) - d
• yj is aligned to a gap, F(i,j) = F(i,j-1) - d
• The best score up to (i,j) will be the smallest of the three
options
Michael
Michael
Schatz
Schatz
Michael Schatz
Smith-Waterman is OPTIMAL but computationally slow
• SW search requires computing of matrix of
scores at every possible alignment position with
every possible gap.
• Compute task increases with the product of the
lengths of two sequence to be compared (N2)
• Difficult for comparison of one small sequence
to a much larger one, very difficult for two large
sequences, essentially impossible to search very
large databases.
Scoring Similarity
1) Can only score aligned sequences
2) DNA is usually scored as identical or not
3) modified scoring for gaps - single vs. multiple base
gaps (gap extension)
4) Protein AAs have varying degrees of similarity
– a. # of mutations to convert one to another
– b. chemical similarity
– c. observed mutation frequencies
5) PAM matrix calculated from observed mutations in
protein families
Protein Scoring Systems
• Amino acids have different biochemical and physical properties
that influence their relative replaceability in evolution.
tiny
aliphatic
P
C S+S
I
V
A
L
hydrophobic M
Y
F
small
G
G
CSH
T
S
D
K
W
H
N
E
R
Q
aromatic
positive
polar
charged
The PAM 250 scoring matrix
PAM Vs. BLOSUM
PAM100
PAM120
PAM160
PAM200
PAM250
=
=
=
=
=
BLOSUM90
BLOSUM80
BLOSUM60
BLOSUM52
BLOSUM45
More distant sequences
PAM120 for general use
PAM60 for close relations
PAM250 for distant relations

BLOSUM62 for general use
BLOSUM80 for close relations
BLOSUM45 for distant relations

Search with Protein, not DNA
Sequences
1) 4 DNA bases vs. 20 amino acids - less
chance similarity
2) can have varying degrees of similarity
between different AAs
- # of mutations, chemical similarity, PAM matrix
3) protein databanks are much smaller than
DNA databanks
FASTA
1) A faster method to find similarities and
make alignments – capable of searching
many sequences (an entire database)
2) Only searches near the diagonal of the
alignment matrix
3) Produces a statistic for each alignment
(more on this in the next lecture)
FASTA
1) Derived from logic of the dot plot
– compute best diagonals from all frames of alignment
2) Word method looks for exact matches between
words in query and test sequence
–
–
–
–
–
hash tables (fast computer technique)
Only matches exactly identical words
DNA words are usually 6 bases
protein words are 1 or 2 amino acids
only searches for diagonals in region of word
matches = faster searching
FASTA Format
• simple format used by almost all programs
• >header line with a [return] at end
• Sequence (no specific requirements for line
length, characters, etc)
>URO1 uro1.seq
Length: 2018
November 9, 2000 11:50
Type: N
Check: 3854
CGCAGAAAGAGGAGGCGCTTGCCTTCAGCTTGTGGGAAATCCCGAAGATGGCCAAAGACA
ACTCAACTGTTCGTTGCTTCCAGGGCCTGCTGATTTTTGGAAATGTGATTATTGGTTGTT
GCGGCATTGCCCTGACTGCGGAGTGCATCTTCTTTGTATCTGACCAACACAGCCTCTACC
CACTGCTTGAAGCCACCGACAACGATGACATCTATGGGGCTGCCTGGATCGGCATATTTG
TGGGCATCTGCCTCTTCTGCCTGTCTGTTCTAGGCATTGTAGGCATCATGAAGTCCAGCA
GGAAAATTCTTCTGGCGTATTTCATTCTGATGTTTATAGTATATGCCTTTGAAGTGGCAT
CTTGTATCACAGCAGCAACACAACAAGACTTTTTCACACCCAACCTCTTCCTGAAGCAGA
TGCTAGAGAGGTACCAAAACAACAGCCCTCCAAACAATGATGACCAGTGGAAAAACAATG
GAGTCACCAAAACCTGGGACAGGCTCATGCTCCAGGACAATTGCTGTGGCGTAAATGGTC
CATCAGACTGGCAAAAATACACATCTGCCTTCCGGACTGAGAATAATGATGCTGACTATC
CCTGGCCTCGTCAATGCTGTGTTATGAACAATCTTAAAGAACCTCTCAACCTGGAGGCTT
..
FASTA Algorithm
Makes Longest Diagonal
3) after all diagonals found, tries to join
diagonals by adding gaps
4) computes alignments in regions of
best diagonals
FASTA Alignments
FASTA Results - Alignment
SCORES
Init1: 1515 Initn: 1565 Opt: 1687 z-score: 1158.1 E(): 2.3e-58
>>GB_IN3:DMU09374
(2038 nt)
initn: 1565 init1: 1515 opt: 1687 Z-score: 1158.1 expect(): 2.3e-58
66.2% identity in 875 nt overlap
(83-957:151-1022)
60
70
80
90
100
110
u39412.gb_pr CCCTTTGTGGCCGCCATGGACAATTCCGGGAAGGAAGCGGAGGCGATGGCGCTGTTGGCC
|| ||| | ||||| |
||| |||||
DMU09374
AGGCGGACATAAATCCTCGACATGGGTGACAACGAACAGAAGGCGCTCCAACTGATGGCC
130
140
150
160
170
180
120
130
140
150
160
170
u39412.gb_pr GAGGCGGAGCGCAAAGTGAAGAACTCGCAGTCCTTCTTCTCTGGCCTCTTTGGAGGCTCA
|||||||||
|| |||
|
| || ||| |
|| || ||||| ||
DMU09374
GAGGCGGAGAAGAAGTTGACCCAGCAGAAGGGCTTTCTGGGATCGCTGTTCGGAGGGTCC
190
200
210
220
230
240
180
190
200
210
220
230
u39412.gb_pr TCCAAAATAGAGGAAGCATGCGAAATCTACGCCAGAGCAGCAAACATGTTCAAAATGGCC
||| | ||||| ||
|||
||||
| || | |||||||| || ||| ||
DMU09374
AACAAGGTGGAGGACGCCATCGAGTGCTACCAGCGGGCGGGCAACATGTTTAAGATGTCC
250
260
270
280
290
300
240
250
260
270
280
290
u39412.gb_pr AAAAACTGGAGTGCTGCTGGAAACGCGTTCTGCCAGGCTGCACAGCTGCACCTGCAGCTC
||||||||||
||||| |
|||||| |||| |||
|| ||| || |
DMU09374
AAAAACTGGACAAAGGCTGGGGAGTGCTTCTGCGAGGCGGCAACTCTACACGCGCGGGCT
310
320
330
340
350
360
Summary
• Identity, similarity, homology
• Analyze sequence similarity by dotplots
• window/stringency
•
•
•
•
•
•
Alignment of text strings by edit distance
Scoring of aligned amino acids
Gap penalties
Global vs. local alignment
Dynamic Programming (Smith Waterman)
FASTA method
Next Lecture: Searching Sequence Databases