SLIDE 1
Modern Computational Statistics Lecture 20: Applications in - - PowerPoint PPT Presentation
Modern Computational Statistics Lecture 20: Applications in - - PowerPoint PPT Presentation
Modern Computational Statistics Lecture 20: Applications in Computational Biology Cheng Zhang School of Mathematical Sciences, Peking University December 09, 2019 Introduction 2/23 While modern statistical approaches have been quite
SLIDE 2
SLIDE 3
Challenges in Computational Biology
3/23 Adapted from Narges Razavian 2013
SLIDE 4
Phylogenetic Inference
4/23
The goal of phylogenetic inference is to reconstruct the evolution history (e.g., phylogenetic trees) from molecular sequence data (e.g., DNA, RNA or protein sequences)
Molecular Sequence Data
Taxa Species A Species B Species C Species D Characters ATGAACAT ATGCACAC ATGCATAT ATGCATGC
Phylogenetic Tree
D A B C
Lots of modern biological and medical applications: predict the evolution of influenza viruses and help vaccine design, etc.
SLIDE 5
Example: B Cell Evolution
5/23
This happens inside of you!
SLIDE 6
Example: B Cell Evolution
5/23
This happens inside of you!
SLIDE 7
Example: B Cell Evolution
5/23
This happens inside of you!
SLIDE 8
Example: B Cell Evolution
5/23
This happens inside of you! These inferences guide rational vaccine design.
SLIDE 9
Bayesian Phylogenetics
6/23
ATGAAC · · · ATGCAC · · · ATGCAT · · · ATGCAT · · ·
(τ, q)
y1y2y3y4y5y6
SLIDE 10
Bayesian Phylogenetics
6/23
pa ch
ATGAAC · · · ATGCAC · · · ATGCAT · · · ATGCAT · · ·
(τ, q)
y1y2y3y4y5y6
e
Evolution model: p(ch|pa, qe) qe: amount of evolution on e.
SLIDE 11
Bayesian Phylogenetics
6/23
A A A A
ATGAAC · · · ATGCAC · · · ATGCAT · · · ATGCAT · · ·
(τ, q)
y1y2y3y4y5y6
Evolution model: p(ch|pa, qe) qe: amount of evolution on e. Likelihood p(Y |τ, q) = η(ai
ρ)
- (u,v)∈E(τ)
Pai
uai v(quv)
SLIDE 12
Bayesian Phylogenetics
6/23
A A A A
ATGAAC · · · ATGCAC · · · ATGCAT · · · ATGCAT · · ·
(τ, q)
y1y2y3y4y5y6
Evolution model: p(ch|pa, qe) qe: amount of evolution on e. Likelihood p(Y |τ, q) =
- ai
η(ai
ρ)
- (u,v)∈E(τ)
Pai
uai v(quv)
SLIDE 13
Bayesian Phylogenetics
6/23
T T T T
ATGAAC · · · ATGCAC · · · ATGCAT · · · ATGCAT · · ·
(τ, q)
y1y2y3y4y5y6
Evolution model: p(ch|pa, qe) qe: amount of evolution on e. Likelihood p(Y |τ, q) =
M
- i=1
- ai
η(ai
ρ)
- (u,v)∈E(τ)
Pai
uai v(quv)
SLIDE 14
Bayesian Phylogenetics
6/23
G G G G
ATGAAC · · · ATGCAC · · · ATGCAT · · · ATGCAT · · ·
(τ, q)
y1y2y3y4y5y6
Evolution model: p(ch|pa, qe) qe: amount of evolution on e. Likelihood p(Y |τ, q) =
M
- i=1
- ai
η(ai
ρ)
- (u,v)∈E(τ)
Pai
uai v(quv)
SLIDE 15
Bayesian Phylogenetics
6/23
A C C C
ATGAAC · · · ATGCAC · · · ATGCAT · · · ATGCAT · · ·
(τ, q)
y1y2y3y4y5y6
Evolution model: p(ch|pa, qe) qe: amount of evolution on e. Likelihood p(Y |τ, q) =
M
- i=1
- ai
η(ai
ρ)
- (u,v)∈E(τ)
Pai
uai v(quv)
SLIDE 16
Bayesian Phylogenetics
6/23
A A A A
ATGAAC · · · ATGCAC · · · ATGCAT · · · ATGCAT · · ·
(τ, q)
y1y2y3y4y5y6
Evolution model: p(ch|pa, qe) qe: amount of evolution on e. Likelihood p(Y |τ, q) =
M
- i=1
- ai
η(ai
ρ)
- (u,v)∈E(τ)
Pai
uai v(quv)
SLIDE 17
Bayesian Phylogenetics
6/23
C C T T
ATGAAC · · · ATGCAC · · · ATGCAT · · · ATGCAT · · ·
(τ, q)
y1y2y3y4y5y6
Evolution model: p(ch|pa, qe) qe: amount of evolution on e. Likelihood p(Y |τ, q) =
M
- i=1
- ai
η(ai
ρ)
- (u,v)∈E(τ)
Pai
uai v(quv)
SLIDE 18
Bayesian Phylogenetics
6/23
C C T T
ATGAAC · · · ATGCAC · · · ATGCAT · · · ATGCAT · · ·
(τ, q)
y1y2y3y4y5y6
Evolution model: p(ch|pa, qe) qe: amount of evolution on e. Likelihood p(Y |τ, q) =
M
- i=1
- ai
η(ai
ρ)
- (u,v)∈E(τ)
Pai
uai v(quv)
SLIDE 19
Bayesian Phylogenetics
6/23
C C T T
ATGAAC · · · ATGCAC · · · ATGCAT · · · ATGCAT · · ·
(τ, q)
y1y2y3y4y5y6
Evolution model: p(ch|pa, qe) qe: amount of evolution on e. Likelihood p(Y |τ, q) =
M
- i=1
- ai
η(ai
ρ)
- (u,v)∈E(τ)
Pai
uai v(quv)
Given a proper prior distribution p(τ, q), the posterior is p(τ, q|Y ) ∝ p(Y |τ, q)p(τ, q).
SLIDE 20
Markov chain Monte Carlo
7/23
Random-walk MCMC (MrBayes, BEAST): ◮ simple random perturbation (e.g., Nearest Neighborhood Interchange) to generate new state. NNI Challenges for MCMC ◮ Large search space: (2n − 5)!! unrooted trees (n taxa) ◮ Intertwined parameter space, low acceptance rate, hard to scale to data sets with many sequences.
SLIDE 21
Variational Inference
8/23
q∗(θ) p(θ|x) Q q∗(θ) = arg min
q∈Q
KL (q(θ)p(θ|x)) ◮ VI turns inference into optimization ◮ Specify a variational family of distributions over the model parameters Q = {qφ(θ); φ ∈ Φ} ◮ Fit the variational parameters φ to minimize the distance (often in terms of KL divergence) to the exact posterior
SLIDE 22
Evidence Lower Bound
9/23
L(θ) = Eq(θ)(log p(x, θ)) − Eq(θ)(log q(θ)) ≤ log p(x) ◮ KL is intractable; maximizing the evidence lower bound (ELBO) instead, which only requires the joint probability p(x, θ).
◮ The ELBO is a lower bound on log p(x). ◮ Maximizing the ELBO is equivalent to minimizing the KL.
◮ The ELBO strikes a balance between two terms
◮ The first term encourages q to focus probability mass where the model puts high probability. ◮ The second term encourages q to be diffuse.
◮ As an optimization approach, VI tends to be faster than MCMC, and is easier to scale to large data sets (via stochastic gradient ascent)
SLIDE 23
Subsplit Bayesian Networks
10/23
Inspired by previous works (H¨
- hna and Drummond 2012,
Larget 2013), we can decompose trees into local structures and encode the tree topology space via Bayesian networks!
D A B C A B C D
SLIDE 24
Subsplit Bayesian Networks
10/23
Inspired by previous works (H¨
- hna and Drummond 2012,
Larget 2013), we can decompose trees into local structures and encode the tree topology space via Bayesian networks!
ABC D AB CD
D A B C A B C D
SLIDE 25
Subsplit Bayesian Networks
10/23
Inspired by previous works (H¨
- hna and Drummond 2012,
Larget 2013), we can decompose trees into local structures and encode the tree topology space via Bayesian networks!
ABC D A BC
D
1.0 AB CD A B C D
D A B C A B C D
SLIDE 26
Subsplit Bayesian Networks
10/23
Inspired by previous works (H¨
- hna and Drummond 2012,
Larget 2013), we can decompose trees into local structures and encode the tree topology space via Bayesian networks!
ABC D A BC
D A
B C
D D
1.0 1.0 1.0 1.0 AB CD A B C D
A B C D
1.0 1.0 1.0 1.0
D A B C A B C D
SLIDE 27
Subsplit Bayesian Networks
10/23
Inspired by previous works (H¨
- hna and Drummond 2012,
Larget 2013), we can decompose trees into local structures and encode the tree topology space via Bayesian networks!
S4 S5 S6 S7 S2 S3 S1
ABC D A BC
D A
B C
D D
1.0 1.0 1.0 1.0 AB CD A B C D
A B C D
1.0 1.0 1.0 1.0
D A B C A B C D
SLIDE 28
Probability Estimation Over Tree Topologies
11/23
S4 S5 S6 S7 S2 S3 S1
ABC D A BC
D A
B C
D D
1.0 1.0 1.0 1.0 AB CD A B C D
A B C D
1.0 1.0 1.0 1.0
D A B C A B C D
Rooted Trees psbn(T = τ) = p(S1 = s1)
- i>1
p(Si = si|Sπi = sπi).
SLIDE 29
Probability Estimation Over Tree Topologies
11/23
A B C D
1 2 4 5 3 A B C D
1
r
- t
/ u n r
- t
A B C D
3
r
- t
/ u n r
- t
A A B
C D
A
B CD A BCD
A B C D
A B C D AB CD
S4 S5 S6 S7 S2 S3 S1
Unrooted Trees: psbn(T u = τ) =
- s1∼τ
p(S1 = s1)
- i>1
p(Si = si|Sπi = sπi).
SLIDE 30
Tree Probability Estimation via SBNs
12/23
SBNs can be used to learn a probability distribution based on a collection of trees T = {T1, · · · , TK}. Tk = {Si = si,k, i ≥ 1}, k = 1, . . . , K
Rooted Trees
◮ Maximum Likelihood Estimates: relative frequencies. ˆ pMLE(S1 = s1) = ms1 K , ˆ pMLE(Si = si|Sπi = ti) = msi,ti
- s∈Ci ms,ti
Unrooted Trees
◮ Expectation Maximization ˆ pEM,(n+1) = arg max
p
Ep(S1|T,ˆ
pEM,(n))
- log p(S1) +
- i>1
log p(Si|Sπi)
SLIDE 31
Example: Phylogenetic Posterior Estimation
13/23
10
8
10
6
10
4
10
2
100
log(ground truth)
10
8
10
7
10
6
10
5
10
4
10
3
10
2
10
1
100
log(estimated probability)
CCD
peak 1 peak 2
10
8
10
6
10
4
10
2
100
log(ground truth)
10
8
10
7
10
6
10
5
10
4
10
3
10
2
10
1
100
log(estimated probability)
SBN-EM
peak 1 peak 2
104 105
number of samples
10
2
10
1
100
KL divergence
DS1
ccd sbn-sa sbn-em sbn-em- srf
[Zhang and Matsen, NeurIPS 2018]
◮ Compared to a previous method CCD (Larget, 2013), SBNs significantly reduce the biases for both high probability and low probability trees. ◮ SBNs perform better in the weak data regime.
SLIDE 32
Example: Phylogenetic Posterior Estimation
14/23
Data set (#Taxa, #Sites) Tree space size Sampled trees KL divergence to ground truth SRF CCD SBN-SA SBN-EM SBN-EM-α DS1 (27, 1949) 5.84×1032 1228 0.0155 0.6027 0.0687 0.0136 0.0130 DS2 (29, 2520) 1.58×1035 7 0.0122 0.0218 0.0218 0.0199 0.0128 DS3 (36, 1812) 4.89×1047 43 0.3539 0.2074 0.1152 0.1243 0.0882 DS4 (41, 1137) 1.01×1057 828 0.5322 0.1952 0.1021 0.0763 0.0637 DS5 (50, 378) 2.84×1074 33752 11.5746 1.3272 0.8952 0.8599 0.8218 DS6 (50, 1133) 2.84×1074 35407 10.0159 0.4526 0.2613 0.3016 0.2786 DS7 (59, 1824) 4.36×1092 1125 1.2765 0.3292 0.2341 0.0483 0.0399 DS8 (64, 1008) 1.04×10103 3067 2.1653 0.4149 0.2212 0.1415 0.1236
[Zhang and Matsen, NeurIPS 2018]
Remark: Unlike previous methods, SBNs are flexible enough to provide accurate approximations to real data posteriors!
SLIDE 33
Variational Bayesian Phylogenetic Inference
15/23
◮ Approximating Distribution:
tree topology
Qφ(τ)
SLIDE 34
Variational Bayesian Phylogenetic Inference
15/23
◮ Approximating Distribution: Qφ,ψ(τ, q)
tree topology
Qφ(τ) ·
branch length
Qψ(q|τ)
SLIDE 35
Variational Bayesian Phylogenetic Inference
15/23
◮ Approximating Distribution: Qφ,ψ(τ, q)
tree topology
Qφ(τ) ·
branch length
Qψ(q|τ) ◮ Multi-sample Lower Bound:
LK(φ, ψ) = EQφ,ψ(τ 1:K,q1:K) log
- 1
K
K
- i=1
p(Y |τ i, qi)p(τ i, qi) Qφ(τ i)Qψ(qi|τ i)
SLIDE 36
Variational Bayesian Phylogenetic Inference
15/23
◮ Approximating Distribution: Qφ,ψ(τ, q)
tree topology
Qφ(τ) ·
branch length
Qψ(q|τ) ◮ Multi-sample Lower Bound:
LK(φ, ψ) = EQφ,ψ(τ 1:K,q1:K) log
- 1
K
K
- i=1
p(Y |τ i, qi)p(τ i, qi) Qφ(τ i)Qψ(qi|τ i)
- ◮ Use stochastic gradient ascent (SGA) to maximize the
lower bound: ˆ φ, ˆ ψ = arg max
φ,ψ
LK(φ, ψ)
◮ φ: VIMCO/RWS ◮ ψ: The Reparameterization Trick
SLIDE 37
Structured Parameterization
16/23
SBNs Parameters
p(S1 = s1) = exp(φs1)
- sr∈Sr exp(φsr),
p(Si = s|Sπi = t) = exp(φs|t)
- s∈S·|t exp(φs|t)
Branch Length Parameters
Qψ(q|τ) =
- e∈E(τ)
pLognormal (qe | µ(e, τ), σ(e, τ)) ◮ Simple Split µs(e, τ) = ψµ
e/τ, σs(e, τ) = ψσ e/τ.
W Z e
ψµ(W, Z)
µs(e, τ)
SLIDE 38
Structured Parameterization
16/23
SBNs Parameters
p(S1 = s1) = exp(φs1)
- sr∈Sr exp(φsr),
p(Si = s|Sπi = t) = exp(φs|t)
- s∈S·|t exp(φs|t)
Branch Length Parameters
Qψ(q|τ) =
- e∈E(τ)
pLognormal (qe | µ(e, τ), σ(e, τ)) ◮ Simple Split µs(e, τ) = ψµ
e/τ, σs(e, τ) = ψσ e/τ.
◮ Primary Subsplit Pair (PSP) µpsp(e, τ) = ψµ
e/τ +
- s∈e/
/τ ψµ s
σpsp(e, τ) = ψσ
e/τ +
- s∈e/
/τ ψσ s .
W
ψµ ( W1 , W2 | W , Z )
Z
ψµ ( Z1 , Z2 | W , Z )
e
ψµ(W, Z)
W1 W2 Z1 Z2
+
µpsp(e, τ)
SLIDE 39
Stochastic Gradient Estimators
17/23
SBNs Parameters φ. With τ j, qj
iid
∼ Qφ,ψ(τ, q) ◮ VIMCO. [Minh and Rezende, ICML 2016] ∇φLK(φ, ψ) ≃
K
- j=1
- ˆ
LK
j|−j(φ, ψ) − ˜
wj ∇φ log Qφ(τ j). ◮ RWS. [Bornschein and Bengio, ICLR 2015] ∇φLK(φ, ψ) ≃
K
- j=1
˜ wj∇φ log Qφ(τ j). Branch Length Parameters ψ. gψ(ǫ|τ) = exp(µψ,τ + σψ,τ ⊙ ǫ). ◮ Reparameterization Trick. Let fφ,ψ(τ, q) = p(Y |τ,q)p(τ,q)
Qφ(τ)Qψ(q|τ).
∇ψLK(φ, ψ) ≃
K
- j=1
˜ wj∇ψ log fφ,ψ(τ j, gψ(ǫj|τ j)) where τ j
iid
∼ Qφ(τ), ǫj
iid
∼ N(0, I).
SLIDE 40
The VBPI Pipeline
18/23
Qφ(τ)
SLIDE 41
The VBPI Pipeline
18/23
Ancestral sampling for SBNs
S4 S5 S6 S7 S2 S3 S1
D A B C A B C D
SLIDE 42
The VBPI Pipeline
18/23
Ancestral sampling for SBNs
S4 S5 S6 S7 S2 S3 S1
ABC D AB CD
D A B C A B C D
SLIDE 43
The VBPI Pipeline
18/23
Ancestral sampling for SBNs
S4 S5 S6 S7 S2 S3 S1
ABC D A BC
D
1.0 AB CD A B C D
D A B C A B C D
SLIDE 44
The VBPI Pipeline
18/23
Ancestral sampling for SBNs
S4 S5 S6 S7 S2 S3 S1
ABC D A BC
D A
B C
D D
1.0 1.0 1.0 1.0 AB CD A B C D
A B C D
1.0 1.0 1.0 1.0
D A B C A B C D
SLIDE 45
The VBPI Pipeline
18/23
Qφ(τ)
sample
e.g., ancestral sampling for SBNs
B D A C
τ 1
C D A B
τ 2
. . .
B C A D
τ K
SLIDE 46
The VBPI Pipeline
18/23
Qφ(τ)
sample
e.g., ancestral sampling for SBNs
B D A C
τ 1
C D A B
τ 2
. . .
B C A D
τ K
Qψ(q|τ)
SLIDE 47
The VBPI Pipeline
18/23
Qφ(τ)
sample
e.g., ancestral sampling for SBNs
B D A C
τ 1
C D A B
τ 2
. . .
B C A D
τ K
Qψ(q|τ)
sample
e.g., Lognormal for branch lengths
B D A C
(τ 1, q1)
C D A B
(τ 2, q2)
. . .
B C A D
(τ K, qK)
SLIDE 48
The VBPI Pipeline
18/23
Qφ(τ)
sample
e.g., ancestral sampling for SBNs
B D A C
τ 1
C D A B
τ 2
. . .
B C A D
τ K
Qψ(q|τ)
sample
e.g., Lognormal for branch lengths
B D A C
(τ 1, q1)
C D A B
(τ 2, q2)
. . .
B C A D
(τ K, qK)
LK(φ, ψ)
multi-sample lower bound
SLIDE 49
The VBPI Pipeline
18/23
Qφ(τ)
sample
e.g., ancestral sampling for SBNs
B D A C
τ 1
C D A B
τ 2
. . .
B C A D
τ K
Qψ(q|τ)
sample
e.g., Lognormal for branch lengths
B D A C
(τ 1, q1)
C D A B
(τ 2, q2)
. . .
B C A D
(τ K, qK)
LK(φ, ψ)
multi-sample lower bound
SGA update
SLIDE 50
Performance on Synthetic Data
19/23
A simulated study on unrooted phylogenetic trees with 8 leaves (10395 trees). The target distribution is a random sample from the symmetric Dirichlet distribution Dir(β1), β = 0.008
50 100 150 200
Thousand Iterations
3.0 2.5 2.0 1.5 1.0 0.5 0.0
Evidence Lower Bound
EXACT VIMCO(20) VIMCO(50) RWS(20) RWS(50)
0.02 0.00
50 100 150 200
Thousand Iterations
10
1
100
KL Divergence
VIMCO(20) VIMCO(50) RWS(20) RWS(50)
10
4
10
3
10
2
10
1
Ground truth
10
4
10
3
10
2
10
1
Variational approximation
VIMCO(50) RWS(50)
[Zhang and Matsen, ICLR 2019]
ELBOs approach 0 quickly ⇒ SBNs approximations are flexible. More samples in the multi-sample ELBOs could be helpful.
SLIDE 51
Performance on Real Data
20/23
50 100 150 200
Thousand Iterations
10
1
100 101
KL Divergence
VIMCO(10) VIMCO(20) RWS(10) RWS(20) MCMC
50 100 150 200
Thousand Iterations
10
1
100 101
KL Divergence
VIMCO(10) + PSP VIMCO(20) + PSP RWS(10) + PSP RWS(20) + PSP MCMC
7042 7040 7038 7036
GSS
7042 7041 7040 7039 7038 7037 7036
VBPI
[Zhang and Matsen, ICLR 2019]
◮ More samples ⇒ better exploration ⇒ better approximation ◮ More flexible branch length distributions across tree topologies (PSP) ease training and improve approximation ◮ Outperform MCMC via much more efficient tree space exploration and branch length updates
SLIDE 52
Performance on Real Data
21/23
Data set Marginal Likelihood (NATs) VIMCO(10) VIMCO(20) VIMCO(10)+PSP VIMCO(20)+PSP SS DS1
- 7108.43(0.26)
- 7108.35(0.21)
- 7108.41(0.16)
- 7108.42(0.10)
- 7108.42(0.18)
DS2
- 26367.70(0.12)
- 26367.71(0.09)
- 26367.72(0.08)
- 26367.70(0.10)
- 26367.57(0.48)
DS3
- 33735.08(0.11)
- 33735.11(0.11)
- 33735.10(0.09)
- 33735.07(0.11)
- 33735.44(0.50)
DS4
- 13329.90(0.31)
- 13329.98(0.20)
- 13329.94(0.18)
- 13329.93(0.22)
- 13330.06(0.54)
DS5
- 8214.36(0.67)
- 8214.74(0.38)
- 8214.61(0.38)
- 8214.55(0.43)
- 8214.51(0.28)
DS6
- 6723.75(0.68)
- 6723.71(0.65)
- 6724.09(0.55)
- 6724.34(0.45)
- 6724.07(0.86)
DS7
- 37332.03(0.43)
- 37331.90(0.49)
- 37331.90(0.32)
- 37332.03(0.23)
- 37332.76(2.42)
DS8
- 8653.34(0.55)
- 8651.54(0.80)
- 8650.63(0.42)
- 8650.55(0.46)
- 8649.88(1.75)
[Zhang and Matsen, ICLR 2019]
◮ Competitive to state-of-the-art (stepping-stone), dramatically reducing cost at test time: VBPI(1000) vs SS(100,000) ◮ PSP alleviates the demand for large samples, reducing computation while maintaining approximation accuracy
SLIDE 53
Conclusion
22/23
◮ We introduced VBPI, a general variational framework for Bayesian phylogenetic inference. ◮ VBPI allows efficient learning on both tree topology and branch lengths, providing competitive performance to MCMC while requiring much less computation. ◮ Can be used for further statistical analysis (e.g., marginal likelihood estimation) via importance sampling. ◮ There are many extensions, including more flexible branch length distributions, more general models, designing adaptive transition kernels in MCMC approaches, etc.
SLIDE 54
References
23/23
◮ Sebastian H´
- hna and Alexei J. Drummond. Guided tree
topology proposals for Bayesian phyloge- netic inference.
- Syst. Biol., 61(1):1–11, January 2012.