new Some sequencing technologies
CS262 Lecture 12, Win07, Batzoglou
Molecular Inversion Probes
CS262 Lecture 12, Win07, Batzoglou
Single Molecule Array for Genotyping--Solexa
CS262 Lecture 12, Win07, Batzoglou
Nanopore Sequencing
CS262 Lecture 12, Win07, Batzoglou
http://www.mcb.harvard.edu/branton/index.htm
Pyrosequencing on a chip
Mostafa Ronaghi, Stanford Genome Technologies Center
CS262 Lecture 12, Win07, Batzoglou
454 Life Sciences
Polony Sequencing
CS262 Lecture 12, Win07, Batzoglou
Some future directions for sequencing
1. Personalized genome sequencing
Find your ~3,000,000 single nucleotide polymorphisms (SNPs) Find your rearrangements Goals:
Link genome with phenotype Provide personalized diet and medicine (???) designer babies, big-brother insurance companies
Timeline:
Inexpensive sequencing: Genotypephenotype association: Personalized drugs: 2010-2015 2010-??? 2015-???
CS262 Lecture 12, Win07, Batzoglou
Some future directions for sequencing
2. Environmental sequencing
Find your flora:
organisms living in your body
External organs: skin, mucous membranes Gut, mouth, etc.
Normal flora: >200 species, >trillions of individuals Floradisease, floranon-optimal health associations Timeline:
Inexpensive research sequencing: Research & associations Personalized sequencing today within next 10 years 2015+
Find diversity of organisms living in different environments
Hard to isolate Assembly of all organisms at once
CS262 Lecture 12, Win07, Batzoglou
Some future directions for sequencing
1. Organism sequencing
Sequence a large fraction of all organisms Deduce ancestors
Reconstruct ancestral genomes Synthesize ancestral genomes Clone--Jurassic park!
Study evolution of function
Find functional elements within a genome How those evolved in different organisms Find how modules/machines composed of many genes evolved
CS262 Lecture 12, Win07, Batzoglou
1
4
Phylogeny Tree Reconstruction
3
2
5
1
4
2
3
5
Phylogenetic Trees
Nodes: species Edges: time of independent evolution Edge length represents evolution time
AKA genetic distance Not necessarily chronological time
CS262 Lecture 12, Win07, Batzoglou
Inferring Phylogenetic Trees
Trees can be inferred by several criteria: Morphology of the organisms
Can lead to mistakes!
Sequence comparison Example: Orc: Elf: Dwarf: Hobbit: Human:
CS262 Lecture 12, Win07, Batzoglou
ACAGTGACGCCCCAAACGT ACAGTGACGCTACAAACGT CCTGTGACGTAACAAACGA CCTGTGACGTAGCAAACGA CCTGTGACGTAGCAAACGA
Modeling Evolution
During infinitesimal time t, there is not enough time for two substitutions to happen on the same nucleotide So we can estimate P(x | y, t), for x, y {A, C, G, T} Then let P(A|A, t) ...... ... ... P(T|A, t) ...... P(A|T, t) ... ... P(T|T, t)
S(t) =
x t y
CS262 Lecture 12, Win07, Batzoglou
Modeling Evolution
Reasonable assumption: multiplicative (implying a stationary Markov process) S(t+t') = S(t)S(t') That is to say, P(x | y, t+t') = z P(x | z, t) P(z | y, t') Jukes-Cantor: constant rate of evolution 1 - 3 For short time , S( ) = I+R =
CS262 Lecture 12, Win07, Batzoglou
A
C
T
G
1 - 3 1 - 3
1 - 3
Modeling Evolution
Jukes-Cantor:
S(t+) = S(t)S() = S(t)(I + R)
For longer times, r(t) s(t) s(t) s(t) s(t) r(t) s(t) s(t) s(t) s(t) r(t) s(t) s(t) s(t) s(t) r(t)
Therefore, (S(t+) S(t))/ = S(t) R At the limit of 0, S'(t) = S(t) R Equivalently, r' = -3r + 3s s' = -s + r Those diff. equations lead to: r(t) = (1 + 3 e-4t) s(t) = (1 e-4t)
S(t) =
Where we can derive: r(t) = (1 + 3 e-4t) s(t) = (1 e-4t)
CS262 Lecture 12, Win07, Batzoglou
Modeling Evolution
Kimura: Transitions: A/G, C/T Transversions: A/T, A/C, G/T, C/G Transitions (rate ) are much more likely than transversions (rate ) r(t) s(t) u(t) s(t) s(t) r(t) s(t) u(t) u(t) s(t) r(t) s(t) s(t) u(t) s(t) r(t)
S(t) =
Where
s(t) = (1 e-4t) u(t) = (1 + e-4t e-2(+)t) r(t) = 1 2s(t) u(t)
CS262 Lecture 12, Win07, Batzoglou
Phylogeny and sequence comparison
Basic principles: Degree of sequence difference is proportional to length of independent sequence evolution Only use positions where alignment is pretty certain avoid areas with (too many) gaps
CS262 Lecture 12, Win07, Batzoglou
Distance between two sequences
Given sequences xi, xj, Define dij = distance between the two sequences One possible definition: dij = fraction f of sites u where xi[u] xj[u] Better model (Jukes-Cantor): f = 3 s(t) = (1 e-4t) e-4t = f log (e-4t) = log (1 4/3 f) -4t = log(1 4/3 f) dij = t = - -1 log(1 4/3 f)
CS262 Lecture 12, Win07, Batzoglou
A simple clustering method for building tree
UPGMA (unweighted pair group method using arithmetic averages) Or the Average Linkage Method Given two disjoint clusters Ci, Proofsequences, Cj of
Ci,Cl dpq + Cj,Cl dpq 1 dkl = dij = {p Ci, q Cj}dpq | + |C |) |C | (|Ci j l
|Ci| |Cj|
Claim that if Ck = Ci Cj, then distance to l|) Ci,Cl dpqcluster Cl is:|) Cj,Cl dpq |Ci|/(|Ci||C another + |Cj|/(|Cj||Cl dil |Ci| + djl |Cj| dkl = |Ci| + |Cj|
CS262 Lecture 12, Win07, Batzoglou
= (|Ci| + |Cj|) |Ci| dil + |Cj| djl = (|Ci| + |Cj|)
Algorithm: Average Linkage
Initialization:
Assign each xi into its own cluster Ci Define one leaf per sequence, height 0
3 5 1 4
Iteration:
Find two clusters Ci, Cj s.t. dij is min Let Ck = Ci Cj Define node connecting Ci, Cj, & place it at height dij/2 Delete Ci, Cj
2
Termination:
When two clusters i, j remain, place root at height dij/2 CS262 Lecture 12, Win07, Batzoglou
1 4 2 3 5
Example
v v w x y z v v w x yz CS262 Lecture 12, Win07, Batzoglou
0 0
w
6 0
x
8 8 0
y
8 8 4 0
z
8 8 4 2 0
v v w xyz
0
w
6 0
xyz
8 8 0
vw xyz vw xyz
0 8 0
w
6 0
x
8 8 0
yz
8 8 4 0 v w x 3
4
2 1 y z
Ultrametric Distances and Molecular Clock
Definition: A distance function d(.,.) is ultrametric if for any three distances dij dik dij, it is true that dij dik = dij The Molecular Clock: The evolutionary distance between species x and y is 2 the Earth time to reach the nearest common ancestor That is, the molecular clock has constant rate in all species
The molecular clock results in ultrametric distances years
1 4 2 3 5
CS262 Lecture 12, Win07, Batzoglou
Ultrametric Distances & Average Linkage
1
4
2
3
5
Average Linkage is guaranteed to reconstruct correctly a binary tree with ultrametric distances Proof: Exercise
CS262 Lecture 12, Win07, Batzoglou
Weakness of Average Linkage
Molecular clock: all species evolve at the same rate (Earth time)
However, certain species (e.g., mouse, rat) evolve much faster Example where UPGMA messes up: Correct tree AL tree
3 2
1
CS262 Lecture 12, Win07, Batzoglou
4
1
4
2
3
Additive Distances
1 8 3 9 12 4 7 5 10 2 6
1,4
13
d
11
Given a tree, a distance measure is additive if the distance between any pair of leaves is the sum of lengths of edges connecting them Given a tree T & additive distances dij, can uniquely reconstruct edge lengths: Find two neighboring leaves i, j, with common parent k Place parent node k at distance dkm = (dim + djm dij) from any node m i, j
CS262 Lecture 12, Win07, Batzoglou
Additive Distances
z x
y For any four leaves x, y, z, w, consider the three sums d(x, y) + d(z, w) d(x, z) + d(y, w) d(x, w) + d(y, z)
w
One of them is smaller than the other two, which are equal d(x, y) + d(z, w) < d(x, z) + d(y, w) = d(x, w) + d(y, z)
CS262 Lecture 12, Win07, Batzoglou
Reconstructing Additive Distances Given T
x
D
v v w x y z
0
y 4
5 3
T
3 4 6
w
10 0
x
17 15 0
y
16 14 9 0
z
16 14 15 14 0 z 7
w
v
If we know T and D, but do not know the length of each leaf, we can reconstruct those lengths
CS262 Lecture 12, Win07, Batzoglou
Reconstructing Additive Distances Given T
x
D
v v w x y z
0
y
T
w
10 0
x
17 15 0
y
16 14 9 0
z
16 14 15 14 0 v z w
CS262 Lecture 12, Win07, Batzoglou
Reconstructing Additive Distances Given T D
v v w x y z a a x y
CS262 Lecture 12, Win07, Batzoglou
w
0
x
y
z
y
x
0
10 17 16 16 15 14 14 0 9 0 15 14 0 z
T
a
w
D1
x
0
y
9 0
z
15 14 0
dax = (dvx + dwx dvw) day = (dvy + dwy dvw) daz = (dvz + dwz dvw)
v
0
11 10 10
z
Reconstructing Additive Distances Given T D1
a a x y z
0
x
0
y
9 0
z
y 4 z 15 14 0
x 5 b 7 3 c 3 a 4 w
11 10 10
T
D2
a a b z
CS262 Lecture 12, Win07, Batzoglou
b
6 0
z
10 10 0
0
D3
a a c
0
c
3 0
6 d(a, c) = 3 d(b, c) = d(a, b) d(a, c) = 3 d(c, z) = d(a, z) d(a, c) = 7 d(b, x) = d(a, x) d(a, b) = 5 d(b, y) = d(a, y) d(a, b) = 4 d(a, w) = d(z, w) d(a, z) = 4 d(a, v) = d(z, v) d(a, z) = 6 Correct!!!
v
Neighbor-Joining
Guaranteed to produce the correct tree if distance is additive May produce a good tree even when distance is not additive
1 Step 1: Finding neighboring leaves Define Dij = (N 2) dij k i dik k j djk 0.4 0.1 0.1 3 0.1
0.4
2 Claim: The above "magic trick" ensures that Dij is minimal iff i, j are neighbors
4
CS262 Lecture 12, Win07, Batzoglou
Algorithm: Neighbor-joining
Initialization:
Define T to be the set of leaf nodes, one per sequence Let L = T
Iteration:
Pick i, j s.t. Dij is minimal Define a new node k, and set dkm = (dim + djm dij) for all m L Add k to T, with edges of lengths dik = (dij + ri rj), djk = dij dik Remove i, j from L; Add k to L
Termination:
When L consists of two nodes, i, j, and the edge between them of length dij
CS262 Lecture 12, Win07, Batzoglou
Parsimony direct method not using distances
One of the most popular methods:
GIVEN multiple alignment FIND tree & history of substitutions explaining alignment
Idea: Find the tree that explains the observed sequences with a minimal number of substitutions Two computational subproblems: 1. Find the parsimony cost of a given tree (easy) 1. Search through all tree topologies (hard)
CS262 Lecture 12, Win07, Batzoglou
Example: Parsimony cost of one column
{A}
Final cost C = 1
{A} {A, B}
A B A A Cost C+=1
A {A}
CS262 Lecture 12, Win07, Batzoglou
B {B}
A {A}
A {A}
Parsimony Scoring
Given a tree, and an alignment column u Label internal nodes to minimize the number of required substitutions Initialization: Set cost C = 0; node k = 2N 1 (last leaf) Iteration: If k is a leaf, set Rk = { xk[u] } // Rk is simply the character of kth species If k is not a leaf, Let i, j be the daughter nodes; Set Rk = Ri Rj if intersection is nonempty Set Rk = Ri Rj, and C += 1, if intersection is empty Termination: Minimal cost of tree for column u, = C
CS262 Lecture 12, Win07, Batzoglou
Example
{B} {A,B} {A} {A} {A} A {A} A {A} A {A} A {A} B {B}
{B} {A,B}
B {B}
A {B}
CS262 {A}
B Lecture 12, Win07, Batzoglou
Traceback to find ancestral nucleotides
Traceback: 1. Choose an arbitrary nucleotide from R2N 1 for the root 1. Having chosen nucleotide r for parent k, If r Ri choose r for daughter i Else, choose arbitrary nucleotide from Ri
Easy to see that this traceback produces some assignment of cost C
CS262 Lecture 12, Win07, Batzoglou
Example
Admissible with Traceback
x
{A, B} A {A} {A, B} A A
B
Still optimal, but inadmissible with Traceback
B
x
B A B
B
x
A {A} B B A {A} {B} {B} A A A A
B
x
B A B
x
A B
x
A
CS262 Lecture 12, Win07, Batzoglou
B
Probabilistic Methods
xroot t1 x1 x2 t2
A more refined measure of evolution along a tree than parsimony P(x1, x2, xroot | t1, t2) = P(xroot) P(x1 | t1, xroot) P(x2 | t2, xroot) If we use Jukes-Cantor, for example, and x1 = xroot = A, x2 = C, t1 = t2 = 1,
CS262 Lecture 12, Win07, Batzoglou
= pA (1 + 3e-4) (1 e-4) = ()3(1 + 3e-4)(1 e-4)
Probabilistic Methods
xroot
xu x2
x1
xN
If we know all internal labels xu,
P(x1, x2, ..., xN, xN+1, ..., x2N-1 | T, t) = P(xroot)
j rootP(xj
| xparent(j), tj, parent(j))
Usually we don't know the internal labels, therefore
P(x1, x2, ..., xN | T, t) =
CS262 Lecture 12, Win07, Batzoglou
x x
N+1
N+2
...
x
2N-1
P(x1, x2, ..., x2N-1 | T, t)
Felsenstein's Likelihood Algorithm
To calculate P(x1, x2, ..., xN | T, t) Initialization:
Set k = 2N 1 Let P(Lk | a) denote the prob. of all the leaves below node k, given that the residue at k is a
Iteration: Compute P(Lk | a) for all a
If k is a leaf node: Set P(Lk | a) = 1(a = xk) If k is not a leaf node: 1. Compute P(Li | b), P(Lj | b) for all b, for daughter nodes i, j 2. Set P(Lk | a) =
b,c P(b
| a, ti) P(Li | b) P(c | a, tj) P(Lj | c)
Termination: Likelihood at this column = P(x1, x2, ..., xN | T, t) =
CS262 Lecture 12, Win07, Batzoglou
P(L
a
2N-1
| a)P(a)
Probabilistic Methods
Given M (ungapped) alignment columns of N sequences, Define likelihood of a tree:
L(T, t) = P(Data | T, t) =
m=1...M
P(x1m, ..., xnm, T, t)
Maximum Likelihood Reconstruction: Given data X = (xij), find a topology T and length vector t that maximize likelihood L(T, t)
CS262 Lecture 12, Win07, Batzoglou
Current popular methods
HUNDREDS of programs available! http://evolution.genetics.washington.edu/phylip/software.html#methods Some recommended programs: Discrete--Parsimony-based
Rec-1-DCM3 http://www.cs.utexas.edu/users/tandy/mp.html Tandy Warnow and colleagues
Probabilistic
SEMPHY http://www.cs.huji.ac.il/labs/compbio/semphy/ Nir Friedman and colleagues
CS262 Lecture 12, Win07, Batzoglou
Multiple Sequence Alignments
CS262 Lecture 12, Win07, Batzoglou
Protein Phylogenies
Proteins evolve by both duplication and species divergence
CS262 Lecture 12, Win07, Batzoglou
Protein Phylogenies Example
CS262 Lecture 12, Win07, Batzoglou
CS262 Lecture 12, Win07, Batzoglou
CS262 Lecture 12, Win07, Batzoglou
Definition
Given N sequences x1, x2,..., xN:
Insert gaps (-) in each sequence xi, such that All sequences have the same length L Score of the global map is maximum
A faint similarity between two sequences becomes significant if present in many Multiple alignments can point to elements that are conserved among a class of and therefore important in the biology of these organisms The patterns of conservation can help us tell function of the element
CS262 Lecture 12, Win07, Batzoglou
Scoring Function: Sum Of Pairs
Definition: Induced pairwise alignment
A pairwise alignment induced by the multiple alignment
Example: x: y: z: Induces: x: ACGCGG-C; y: ACGC-GAC; x: AC-GCGG-C; z: GCCGC-GAG; y: AC-GCGAG z: GCCGCGAG AC-GCGG-C AC-GC-GAG GCCGC-GAG
CS262 Lecture 12, Win07, Batzoglou
Sum Of Pairs (cont'd)
Heuristic way to incorporate evolution tree:
Hum a n Mo us e Duc k C h ic ke n
Weighted SOP:
CS262 Lecture 12, Win07, Batzoglou
A Profile Representation
T C C C A C G T .6 1 .2 .2 .2 1 .8 .6 A A A A A 1 1 G G G G G G C C C C C T T T T T A A A A A 1 .4 1 .2 T C C T T C C C C C A A A A G .8 .6 .2 .4 .2 .4 .8 .4 1 C C C C T G G G G G G G
Given a multiple alignment M = m1...mn
Replace each column mi with profile entry pi
Frequency of each letter in # gaps Optional: # gap openings, extensions, closings
Can think of this as a "likelihood" of each letter in each position
CS262 Lecture 12, Win07, Batzoglou
Multiple Sequence Alignments
Algorithms
CS262 Lecture 12, Win07, Batzoglou
Multidimensional DP
Generalization of Needleman-Wunsh:
CS262 Lecture 12, Win07, Batzoglou
Multidimensional DP
Example: in 3D (three sequences): 7 neighbors/cell
F(i,j,k) = max{ F(i 1, j 1, k 1) + S(xi, xj, xk), F(i 1, j 1, k ) + S(xi, xj, - ), F(i 1, j , k 1) + S(xi, -, xk), F(i 1, j ,k ) + S(xi, -, - ), F(i , j 1, k 1) + S( -, xj, xk), F(i , j 1, k ) + S( -, xj, - ), F(i , j , k 1) + S( -, -, xk) }
CS262 Lecture 12, Win07, Batzoglou
Multidimensional DP
Running Time: 1. Size of matrix: LN; Where L = length of each sequence N = number of sequences 1. Neighbors/cell: 2N 1 Therefore.............................. O(2N LN)
CS262 Lecture 12, Win07, Batzoglou
Multidimensional DP
How Running Time: do gap states generalize? VERY badly! 1. Size of matrix: LN; Where
Require 2N 1 states, one per combination of gapped/ungapped sequences L = lengthtime: O(2N sequence= O(4N LN) Running of each 2N LN)
Y YZ
N = number of sequences 1. Neighbors/cell: 2N 1
XY
XYZ
Z
Therefore.............................. O(2N LN)
X
CS262 Lecture 12, Win07, Batzoglou
XZ
Progressive Alignment
pxy pxyzw
x y z w
pzw
When evolutionary tree is known: Align closest first, in the order of the tree In each step, align two sequences x, y, or profiles px, py, to generate a new alignment with associated profile presult Weighted version: Tree edges have weights, proportional to the divergence in that edge New profile is a weighted average of two old profiles
CS262 Lecture 12, Win07, Batzoglou
Progressive Alignment
Example
x Profile: (A, C, G, T, -) x y
p = (0.8, 0.2, 0, 0, 0)
y
z
p = (0.6, 0, 0, 0, 0.4)
w
x
y
When evolutionary tree is known:
Align closest first, in the order of the tree xy In each step, align two sequences x, y, or profiles px, py, to generate a new alignment with associated profile presult Result: p = (0...