Robust and Scalable Models of Microbiome Dynamics Travis E. Gibson 1 - - PowerPoint PPT Presentation

robust and scalable models of microbiome dynamics
SMART_READER_LITE
LIVE PREVIEW

Robust and Scalable Models of Microbiome Dynamics Travis E. Gibson 1 - - PowerPoint PPT Presentation

Robust and Scalable Models of Microbiome Dynamics Travis E. Gibson 1 Georg K. Gerber 1 , 2 tgibson@mit.edu ggerber@bwh.harvard.edu 1 Massachusetts Host Microbiome Center Brigham and Womens Hospital & Harvard Medical School 2 Harvard-MIT


slide-1
SLIDE 1

Robust and Scalable Models of Microbiome Dynamics

Travis E. Gibson1

tgibson@mit.edu

Georg K. Gerber1,2

ggerber@bwh.harvard.edu

1Massachusetts Host Microbiome Center

Brigham and Women’s Hospital & Harvard Medical School

2Harvard-MIT Health Sciences and Technology

35th International Conference on Machine Learning July 12th, 2018

Poster #2

slide-2
SLIDE 2

Outline

1 What is the microbiome 2 From sequences to reads to microbial interactions 3 Bayesian nonparametric model for microbe dynamics

2 / 18

slide-3
SLIDE 3

The Microbiome

1 The microbiome is the aggregate of

microorganisms that resides on or within any of a number of human tissues and biofluids:

  • skin, mammary glands, placenta, seminal

fluid, uterus, ovarian follicles, lung, saliva, oral mucosa, conjunctiva, biliary and gastrointestinal tracts) [wikipedia]

2 1014 Microbes in/on your body [Sender et al.

PLoS Biology 2016]

3 3.3 million genes compared to 23,000 human

genes [Qin et al. Nature 2010]

4 Play a role in a variety of human diseases:

  • infections, arthritis, food allergy, cancer,

inflammatory bowel disease, neurological diseases, and obesity/diabetes

3 / 18

slide-4
SLIDE 4

Bacteriotherapy

Bacteriotherapy: communities of bacteria administered to patients for specific therapeutic applications

  • “bugs-as-drugs”

Clostridium difficile infection

  • Causes serious diarrhea (14K deaths/yr)
  • Antibiotics disrupt helpful bacteria in gut
  • Increasingly difficult to treat with conventional therapies (more antibiotics): 20-30%

recurrence rate Pharmacology meets Ecology

  • C. diff

microbial interaction network positive microbe A produces a small molecule (metabolite) that microbe B needs negative two microbes competing for the same niche what if there were 300 bugs in the network?

4 / 18

slide-5
SLIDE 5

Workflow in our lab

batch experiments chemostat animal experiments

  • 16S rRNA on MiSeq (reads)

for relative abundances of species

relative abundance t i m e

  • 16S rRNA qPCR (universal

primers) for bacterial biomass

time total biomass

5 / 18

slide-6
SLIDE 6

Sequencing to obtain microbe relative abundances

  • 16S rRNA on MiSeq (reads)

for relative abundances of species

relative abundance t i m e

  • 16S rRNA qPCR (universal

primers) for bacterial biomass

time total biomass

  • measurements - irregular,

sparse & noisy

@sequence-id TCGCACTCAACGCCCTGCATATGACAAGACAGAATC + <>;##=><9=AAAAAAAAAA9#:<#<;<<<????#=

  • 1. Fastq file
  • 2. Sequences are clustered
  • 3. Read Table

sample1 sample2 ... 11 1004 bacteria1 bacteria2 bacteria3 7 301 275 . . . nucleobase sequence quality scores

Reads ∼ Negative Binomial

6 / 18

slide-7
SLIDE 7

Quantitative PCR for total bacterial biomass

  • 16S rRNA on MiSeq (reads)

for relative abundances of species

relative abundance t i m e

  • 16S rRNA qPCR (universal

primers) for bacterial biomass

time total biomass

  • measurements - irregular,

sparse & noisy

cycle relative fluorescence control templates library

Dilute control template and library for quantification Amplification Plot

1 µg 0.1 µg 10 ng 1 ng 0.1 ng 10 pg 1 pg

initial quantity Ct

Standard Curve

7 / 18

slide-8
SLIDE 8

Irregular and sparse measurements

  • 16S rRNA on MiSeq (reads)

for relative abundances of species

relative abundance t i m e

  • 16S rRNA qPCR (universal

primers) for bacterial biomass

time total biomass

  • measurements - irregular,

sparse & noisy

perturbation time abundance measurements

8 / 18

slide-9
SLIDE 9

Learning microbial interaction networks

  • 16S rRNA on MiSeq (reads)

for relative abundances of species

relative abundance t i m e

  • 16S rRNA qPCR (universal

primers) for bacterial biomass

time total biomass

  • measurements - irregular,

sparse & noisy

Interaction Network Abundance of microbe i at time t : xt,i dxt,i dt = αixt,i + βiix2

t,i +

  • j=i

βijxt,ixt,j growth, self limiting, interaction Previous literature specific to the microbiome

  • Log transform dynamics → Linear

Regression + L2 [Stein et al. PLoS Comput

Biology 2013]

  • Sparse linear regression with bootstrap

aggregation [Fisher et al. PLoS One 2014]

  • Bayesian model with deterministic

dynamics (independent filtering) [Bucci et

  • al. Genome Biology 2016]
  • Extended Kalman Filter [Alshawaqfeh et al.

BMC Genomics 2017]

9 / 18

slide-10
SLIDE 10

Goal with our model and short literature review

Interaction Network

  • 300 species
  • 90,000 interactions
  • Redundant gene

function Interaction Modules Three main contributions in our model

1 Clustering (interaction modules)

  • Dirichlet Process (DP)

[Rasmussen, Advances in Information Processing Systems 2000] [Neal, Journal of computational and graphical statistics, 2000]

2 Edge selection (structure learning, variable

selection)

  • Bayesian Networks

[George and McCulloch, Journal of the ASA, 1993] [Heckerman, A Tutorial on Learning with Bayesian Networks, 2008]

3 Introduction of an auxiliary variable between

the measurement model and latent state

10 / 18

slide-11
SLIDE 11

Back to the basic time series model

  • Abundance of microbe i at time t : xt,i

dxt,i dt = αixt,i + βiix2

t,i +

  • j=i

βijxt,ixt,j + dwt,i dt growth, self limiting, interaction, stochastic disturbance

  • Convert to discrete time

xk+1,i = xk,i +

  • αixk,i + βiix2

k,i +

  • j=i

βijxk,ixk,j

  • ∆k + (wk+1,i − wk,i)
  • ∆k

discrete time step size Next we discuss the three main ingredients to our model

1 Clustering (interaction modules) 2 Edge selection (structure learning, variable selection) 3 Introduction of an auxiliary variable between the measurement model

11 / 18

slide-12
SLIDE 12

Complete Model

Dirichlet Process πc | α ∼ Stick(α) ci | πc ∼ Multinomial(πc) bci,cj | σb ∼ Normal(0, σ2

b)

Edge Selection (Structure) zci,cj | πz ∼ Bernouli(πz) Self Interactions ai,1, ai,2 | σa ∼ Normal(0, σ2

a)

Dynamics xk+1,i | xk, ai, b, c, z, σw ∼ Normal

  • xk,i+xk,i
  • ai,1+ai,2xk,i+

cj=ci

bci,cjzci,cjxk,j

  • , ∆kσ2

w

  • Constraint and Measurement Model

aux qk,i | xk,i ∼ Normal(xk,i, σ2

q)

reads yk,i | qk,i ∼ NegBin(φ(qk), ǫ(qk)) qPCR Qk | qk,i ∼ Normal

iqk,i, σ2 Qk

  • xk,i

qk,i yk,i Qk bℓ,m σb zℓ,m πz ai k ∈ [m] i ∈ [n] ci πc α σa ℓ ∈ Z+ m ∈ Z+ i ∈ [n]

12 / 18

slide-13
SLIDE 13

Simplified model unraveled in time - auxiliary variable

xt+1,i | xt, a ∼ Normal(aT

i f(xt), σ2 xi)

qk,i | xk,i ∼ Normal(xk,i, σ2

q)

qk,i ∼ Uniform[0, L) yk,i | σy, qk,i ∼ Normal≥0(qk,i, σ2

y )

ai ∼ Normal(0, σ2

ai)

x1 x2 x3 · · · xn a q1 q2 q3 · · · qn y1 y2 y3 · · · yn

Prior on q is positive, relaxing the distribution

  • n the dynamics for x

Parameter inference Gibbs step: a(g+1) ∼ pa|x(· | x(g))

  • Direct sampling from the posterior possible (Bayesian Regression!)

Sampling for other variables

  • Collapsed Gibbs sampling for Dirichlet Process and Edge Selection (integrate out a)
  • Filtering is still challenging but easy to design proposals for (MH)

x1 x2 x3 · · · xn a q1 q2 q3 · · · qn y1 y2 y3 · · · yn

13 / 18

slide-14
SLIDE 14

Synthetic experiment

  • Comparing inference with and without clustering enabled
  • Ground truth interaction network

1, 5, 7 9, 11 2, 4, 6, 8 10, 12 3, 13 2 −4 3 −1 1/(abundance time) 1 5 7 9 11 2 4 6 8 10 12 3 13 1 5 7 9 11 2 4 6 8 10 12 3 13

Microbe Interactions (Truth)

2 2 2 2 2 2 2 2 2 2 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3

  • 1
  • 1
  • 1
  • 1
  • 1
  • 1
  • 1
  • 1
  • 1
  • 1
  • 5
  • 5
  • 5
  • 5
  • 5
  • 5
  • 4
  • 4
  • 5
  • 4
  • 4
  • 5
  • 4
  • 4
  • 5
  • 4
  • 4
  • 5
  • 4
  • 4
  • 5
  • 4
  • 4
  • 5
  • 5
  • 5

5

time

1 2 3 4 5

Biological Replicates

100 105 1010

RMSE (log scaling) Forecast Trajectories

Module Learning Off Module Learning On

1 2 3 4 5

Biological Replicates

10 15 20 25

RMSE Interaction Coefficients

C D

14 / 18

slide-15
SLIDE 15

Synthetic experiment continued

1 5 7 9 11 2 4 6 8 10 12 3 13 1 5 7 9 11 2 4 6 8 10 12 3 13

Microbe Co-cluster Proportions

0.2 0.4 0.6 0.8 1

1 5 7 9 11 2 4 6 8 10 12 3 13 1 5 7 9 11 2 4 6 8 10 12 3 13

Microbe Interactions (RMSE=9.49)

  • 5

5

1/(abundance time) 50 100

time

5 10 15

microbiota abundances Forecasted Trajectories (RMSE=1.88)

1 5 7 9 11 2 4 6 8 10 12 3 13 1 5 7 9 11 2 4 6 8 10 12 3 13

Microbe Co-cluster Proportions

1 1 1 1 1 1 1 1 1 1 1 1 1 0.2 0.4 0.6 0.8 1 1 5 7 9 11 2 4 6 8 10 12 3 13 1 5 7 9 11 2 4 6 8 10 12 3 13

Microbe Interactions (RMSE=15.9)

  • 5

5

1/(abundance time) 50 100

time

5 10 15

microbiota abundances Forecasted Trajectories (RMSE=2.06)

A B

15 / 18

slide-16
SLIDE 16

In vivo infection experiment

  • Data from C. difficile challenge experiment performed with 5 germ free mice [Bucci et al.

Genome Biology 2016]

Day 1 Day 28 Day 56

  • C. difficile

GnotoComplex A

13 samples 13 samples

  • C. scindens
  • B. ovatus
  • P. distasonis
  • A. muciniphila
  • R. hominis
  • C. difficile

+ rest

B

−2.1 −0.13 −0.05 1 2 3 4 5 6 7 8 9 10 11 12 13

Microbe Co-cluster Proportions

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 Klebsiella oxytoca 13 Ruminococcus obeum 12 Bacteroides vulgatus 11 Bacteroides fragilis 10 Escherichia coli 9 Proteus mirabilis 8 Clostridium ramosum 7 Clostridium difficile 6 Roseburia hominis 5 Akkermansia muciniphila 4 Parabacteroides distasonis 3 Bacteroides ovatus 2 Clostridium scindens 1 1 2 3 4 5 6 7 8 9 10 11 12 13

Microbe Interaction Strength

  • 14
  • 13
  • 12
  • 11
  • 10
  • 9
  • 8

Klebsiella oxytoca 13 Ruminococcus obeum 12 Bacteroides vulgatus 11 Bacteroides fragilis 10 Escherichia coli 9 Proteus mirabilis 8 Clostridium ramosum 7 Clostridium difficile 6 Roseburia hominis 5 Akkermansia muciniphila 4 Parabacteroides distasonis 3 Bacteroides ovatus 2 Clostridium scindens 1 -

  • +
  • +

+ + + + +

  • +

+ + + + + + + + + + +

  • +

+

  • +

+ + + + + + +

  • +

+ + + + + + + + + + + +

  • +

+ + + + + + + + + + + +

  • +

+

  • +

+ + + + + + + + + +

  • gram/(CFU day)

log10

C

16 / 18

slide-17
SLIDE 17

Three ongoing projects using the model

Marika Ziesack

Silver Lab, Harvard

  • E. coli
  • B. fragilis
  • S. typhimurium
  • B. theta
  • Microbes engineered to overproduces one amino acid
  • Microbes engineered to need three amino acids
  • Compare inference on WT and engineered strains to

prove that engineering was performed. Bryan Hsu

Silver Lab, Harvard

uninfected bacteria bacteriophages infected bacteria

  • Bacteriophages are “bacteria viruses”
  • Using phages as control knobs to modulate microbiome
  • Unlike antibiotics, phages are host specific and can be

present in equilibrium with host

Massachusetts Host Microbiome Center

Transplant

  • Transplant complex bacteria into germ free mice (300+

bacteria)

17 / 18

slide-18
SLIDE 18

Conclusions

1 Dynamical systems model for microbial dynamics based on what we term interaction

modules (probabilistic clusters of latent variables with redundant interaction structure)

2 Bayesian formulation of the stochastic dynamical systems model that propagates

measurement and latent state uncertainty throughout the model

3 Introduction of a temporally varying auxiliary variable technique to enable efficient

inference by relaxing the hard non-negativity constraint on states Time series metagenomics is a growing area for ML applications

  • Metagenomics is a “Big Data” problem, but number of time series samples is still sparse

in comparison to the number of latent dynamical coefficients we need to learn.

Poster #2

18 / 18

slide-19
SLIDE 19

Backup Slides

19 / 18

slide-20
SLIDE 20

Simple example without the intermediate auxiliary variable

xt+1,i | xt, a ∼ Normal≥0(aiTf(xt), σ2

xi)

yt,i | xt,i ∼ Normal≥0(xt,i, σ2

yi)

ai ∼ Normal(0, σ2

ai) x1 x2 x3 · · · xn a y1 y2 y3 · · · yn

Note the truncated dis- tributions for x and y Parameter inference Gibbs step: a(g+1) ∼ pa|x(· | x(g)) pa|x ∝ px|a

Normal≥0(x; µ(a, x), σ2)

pa

Normal(a; 0, σ2)

= e−

1 2σ2 (x−µ(a,x))2

σ √ 2π

  • Φ(∞) − Φ
  • −µ(a, x)

σ e−

1 2σ2 a2

σ √ 2π Sampling for other variables

  • Filtering (sampling from posterior of x) is challenging
  • Can not use collapsed Gibbs sampling for Dirichlet Process or Edge Selection

20 / 18

slide-21
SLIDE 21

Negative Binomial

yk,i | qk ∼ NegBin(φ(qk, rk), ǫ(qk, a0, a1)) φ(qk, rk) = rk qk,i

  • i qk,i

(1) ǫ(qk, a0, a1) = a0 qk,i/

i qk,i

+ a1 (2) NegBin(y; φ, ǫ) =Γ(r + y) y! Γ(r)

  • φ

r + φ y r r + φ r r =1 ǫ

21 / 18

slide-22
SLIDE 22

Priors

  • We use conjugate priors on many variables (e.g., the variance terms (σ2

a, σ2 b, σ2 w) have

Inverse-Chi-squared priors).

  • The module assignments, c, are also updated by a standard Gibbs sampling approach for

Dirichlet Processes.

  • For the concentration parameter α, which has a Gamma prior on α, we use the sampling

method described by Escobar and West 1995.

22 / 18

slide-23
SLIDE 23

Backup: Alternative View of Bacteriotherapy

Interpretable inference for bacteriotherapy design Predict microbial abundances (phenotype) in the presence of the bacteriotherapy

abundance time

Predicted Abundances with new species

bacteriotherapy administered

23 / 18