A generative nonparametric Bayesian model for whole genomes
Authors: Alan Nawzad Amin
Eli N Weinstein
Debora Susan Marks
Venue: NeurIPS 2021 Poster (2021)
Keywords: Bayesian nonparametrics
genomics
hypothesis testing

Abstract: Generative probabilistic modeling of biological sequences has widespread existing and potential use across biology and biomedicine, particularly given advances in high-throughput sequencing, synthesis and editing. However, we still lack methods with nucleotide resolution that are tractable at the scale of whole genomes and that can achieve high predictive accuracy in theory and practice. In this article we propose a new generative sequence model, the Bayesian embedded autoregressive (BEAR) model, which uses a parametric autoregressive model to specify a conjugate prior over a nonparametric Bayesian Markov model. We explore, theoretically and empirically, applications of BEAR models to a variety of statistical problems including density estimation, robust parameter estimation, goodness-of-fit tests, and two-sample tests. We prove rigorous asymptotic consistency results including nonparametric posterior concentration rates. We scale inference in BEAR models to datasets containing tens of billions of nucleotides. On genomic, transcriptomic, and metagenomic sequence data we show that BEAR models provide large increases in predictive performance as compared to parametric autoregressive models, among other results. BEAR models offer a flexible and scalable framework, with theoretical guarantees, for building and critiquing generative models at the whole genome scale.

A Generative Nonparametric Bayesian Model for Whole Genomes
Alan N. Amin' Eli N. Weinstein and Debora S. Marks Program in Systems, Synthetic and Quantitative Biology DDepartment of Systems Biology
Harvard Medical School
Program in Biophysics. Harvard University
Broad Institule of Harvard and MIT
alanamin@fas hww.arr...... edu, eweinsteingg, harvard. debbie@hms. www.aar..d....... edu
Abstract
Generative probabilistic modeling of biological sequences has widespread existing and potential use across biology and biomedicine, particularly given advances in high-throuhhput sequencing, synthesis and editing. However, we still lack methods with nucleotide resolution that are tractable at the scale of whole genomes and that achieve high predictive accuracy in theory and practice. In this article we propose a new generative sequence model, the Bayesian embedded autoregressive (BEAR) model, which uses a parametric autoregressive model to specify a conjugate prior over a nonpara- metric Bayesian Markov model. We explore, theoretically and empirically, applications of BEAR models to a variety of statistical problems includ- ing density estimation, robust parameter estimation, goodness-offli tests, and two-sample tests. We prove rigorous asymptotic consistency results including nonparametrie posterior concentration rates. We scale inference in BEAR models to datasets containing tens of billions of nucleotides. On genomic, transcriptomic, and metagenomic sequence data we show that BEAR models provide large increases in predictive performance as compared parametric autoregressive models, among other results. BEAR models offer flexible and scalable framework, with theoretical guarantees, for building and critiquing generative models at the whole genome scale.
Introduction
Measuring and making DNA is central to modern biology and biomedicine. Generative probabilistic modeling offers a framework for learning from sequencing data and forming experimentally testable predictions of unobserved or future sequences that can be synthesized Existing approaches to genome modcling typically preprocess the data to build a matrix of genetic variants such single nucleotide polymorphisms 25 57
in the laboratory 19 31 03]
However, most modes sequence variation are more complex. Structural variation occurs widely within individuals (e.g. in cancer), between individuals (e.g. in domesticated plant populations) and between species (e.g. in the human microbiome). and methods for detecting and classifying structural variants are heuristic and designed only for predefined types of sequence variation such as repeats [2 [ bi 39 778. Ideally, we would be able to directly model genome sequencing data and/or or assembled genome sequences. However, building gencrative models that work with raw nuclootides, not matrices of alleles, raises the
Thesse authors contributed equally.
35th Conference on Neural Information Processing Systems (NeurIPS 2021).

---

A
B
10(0)
Deakiel
...GTATCC
C
G
Figure I: Overview of the BEAR model. (A) BEAR models employ a Dirichlet prior on Markov transition probabilities that is centered at the prediction of an AR model. (B) De Bruijn graphs showing BEAR transitions with non-zero probability under an example eeta-generating distribution. As the lag L increases, the model has higher resolution.
extreme statistical challenges of having enough flemiblility to account for genomic complexity, interpretability to reach scientific conclusions, and scalability to train billions of nucleotides. Given the relevance genetic analysis to human health, models should also possess strong theoretical guorantees.
Autoregressive (AR) models are a natural starting point for generative genome modcling, since they (1) have been successfully applied to biological sequences, as well as many other types of non- non-biological sequential data, (2) can be designed have interpretable parameters, and (3) can be scaled to big datasets with very long sequences 66 7 However, since AR models are parametric models, they will in general suffer from misspecificatinns as we show empirically in Section G for genomic datasets misspecification can be a serious practical limitation not only for simple AR models but even for deep neural networks.
As an alternative strategy for building generative probabilistic models at the genome scale, we propose in Section the nonparametric "Bayesian embedded autoregressive" (BEAR) model. BEAR models are Bayesian Markov models, with prior on the lag and conjugate Dirichlet priors on the transition probabilities. The hyperparameters of the Dirichlet prior are controlled by an "embedded" AR model with parameters 0 and an overall concentration hyperparameter h, both of which can be optimized via empirical Bayes. Section we show that BEAR models can capture arbitrary data-generating distributions, and establish asymptotic consistency guarantees and convergence rates for nonparametrie density estima- tion. In Section g we show that the optimal h provides diagnostic for whether or not the embedded AR model is misspecified and so by how much, alerting the practitioner when the parameter estimates 0 are untrustworthy. Besides estimation problems, BEAR models can also be used to construct goodness-of-fit tests and two-sample tests, thanks to their analytic marginal likelihoods, and prove consistency results for these tests in Section Finally we apply BEAR models at large scale, to genomic datasets with tens of billions of mucleotides, including whole genome, whole transcriptome, and metagenomic sequencing data; we find that BEAR models can have greatly improved performance over AR models (Section
Crucial to our theoretical and empirical analysis the statistical setting: we assume that the data XI...., Xx consists of finite but possibly variable length strings (with small alphabets) drawn i.i.d. from some underlying distribution p", and study the behavior of estimators and tests N + cc. This setup differs from common theoretical analyses of sequence models outside of biology, which typically consider the limit as the length of an individual sequence goes to infinity [20 In biology, however, observe finite sequences recorded from many individual species, organisms, cells, molecules, etc. and want to generalize to unseen sequences, making N xo the appropriate large data limit.
2 Bayesian embedded autoregressive models
We first. briefly review autoregressive (AR) models applied to sequences of discrete characters. Let /()) denote an autoregressive function with parameter 0 and let. L denote the lag of the autoregressive model; then the AR model generates data as
X, Categorical (0))
(1)

---

where i indexes position in the sequence X and Xi-Li-1 consists of the previous L letters in the sequence. Since sequence length as well as nucleotide or amino acid content is relevant to biological applications, we use start symbol 0 at the beginning and a stop symbol $ at the end of each sequence; letters X; are sampled sequentially starting from the start symbol and continuing until a stop symbol is drawn.
We propose the Bayesian embedded autoregressive (BEAR) model, a Bayesian Markov model that embeds an AR model into its prior. The BEAR model takes the form,
1 fx(()) for all k, h
L ~ ((), Uk Dirichlet X Categorical(xxx....
where (l) is prior on the lag with support up to infinity, h > 0 is concentration hyperparameterr and & is a length L kmer. The BEAR model has three key properties (Fig- D. First, the unrestricted transition parameter v and lag L allow the model to capture exact conditional distributions of p" to arbitrarily high order: p" (X; at L = 1, then (", Xi-1) at L 2, etc.. This property allows the BEAR model to be used for nonparametric density estimation (Section 3. Second, in the limit where h + 0, the BEAR model reduces to the embedded AR model (Eqn. D The optimal h provides measurement of the amount of misspecification in the AR model (Section 4 Third, the choice of the h, 0) to be computed analytically, and (since L is one-dimensional) the total marginal likelihood (((,)) /h, 0) to be estimated tractably. This allows BEAR models to be used for hypothesis testing (Section 5
conjugate Dirichlet prior allows the conditional marginals p((Xn) )
There a variety of ways of performing inference in BEAR models, but for most applications we will focus on empirical Bayes methods that optimize point estimates of L, h and 0. Let #(,bb) denote the number of times the length L kmer k is seen followed by the letter or stop symbol b in the dataset (X.) Using a nucleotide data, KMC, we can compute the count matrix #(,-) for all observed kmers A in tera scraae datasets, even when the matrix does not fit in main memory (Section J.2 [39
kmer counter optimized for
To optimize h and 0, we take advantage of the fact that the log conditional marginal likelhood can be written as a sum over observed kmers,
T (m)) IL (0) + #(k, b)) III. (( TO f + #(k
log (og h,H,0)
Σ log k
(3)
This decomposition lets us construct unbiased stochastic estimates of the gradient with respect to h and 0 by subsampling rows of the count matrix (Section J.1 Empirical Bayes in the BEAR model therefore costs little extra time as compared to standard stochastic gradient-based optimization of the original AR model. Code is available at /ttpp:///gtthub com//ebbbemarrelabbbBee
Toy example We next briefly illustrate the properties and advantages of the BEAR model in simulation. We generated samples from an AR model in which f& (0) depends on & lincarly as a function of both individual positions and pairwise interactions between positions, with the strength of the pairwise interaction weighted by a parameter 3 (Section L.I. We first fit (using maximum likelihood) linear AR model that lacks pairwise terms and is thus misspecified when 3 > O. Since the AR model is misspecified, it does not asymptotically approach the true data-generating distribution p" (Fig. 2 gray). We next computed the posterior of a vanilla BEAR model without the embedded AR in its prior, instead using the Jeffreys prior Uy ad Dirichlet (1/2.... 1/2). approaches the true data generating distribution, since it is a nonparametric model; however, it underperforms the AR model in the low data regime (Fig. 2A, black). Finally, we fit BEAR model with the misspecified lincar AR model embedded, using our empirical Bayes procedure. The BEAR model performs just well as its embedded AR model in the low data regime, just as well as the vanilla model in the high data regime, and better than both at intermediate values (Fig. A blue and yellow).
The vanilla BEAR model asymptotically
When the AR model is well-specifid,, the empirical Bayes estimates of the parameters under the BEAR model match the maximum likelihood estimates of 0 under the AR model
3

---

A. Density estimation
В. Diagnosing misspecification
8 0.00 Well-specified
8- = 0.6
=- 1.0
Low misspecification
High misspecification
BEAR model
KL divergencc
Model
Medel REAR
BEAR Van
log(h)
10 212 2a2 Dataset size
112 "5y2""ooo Dataset size
12 22220 Mooo
10000
Dataset size
Dataset size
Figure 2: BEAR models detect and avoid misspecification without sacrificing small datasct performance. (A) Estimated KL divergence between simulated data- generating distribution p" and model posterior predictive distribution, as function of dataset size N. Five independent simulations were run; thin lines correspond to individual simulations, thick lines to the average across sinulations. (B) The h misspecification diagnostic as a function of dataset size, for varying 3 Dataset sizes at which is close to convergence for 8 = 0.6 (right) and 8 = 1.0 (left) are marked with vertical lines.
nearly exactly (Fig. S7
When the AR model is misspecified, however, the BEAR model provides a warning: the empirical Bayes estimate of h converges to a non-zero value, rather than zero (Fig. 23). This warning emerges early: converges well before the vanilla model starts outperforming the misspecified AR model.
Related Work The key idea behind BEAR models is to nonparametrically perturb a 48 following similar strategy to the Polya tree method proposed by Berger and Guglielmi [ As in Berger and Guglielmi G we use Dirichlet, priors centered at the parametric model's predictions, and construct tractable goodness-of-ft tests by exploiting Dirichlet-categorical conjugacy. BEAR models extend these ideas from one-dimensional continuous data to finite-lnngth sequences of discrete characters.
parametric model
Markov and AR models have a long history and wide range of applications in biological Compression methods, in particular, often rely accurate density estimation and use Markov or AR models to achieve it [17
sequence analysis 21l 51 5
b
5i
Girl
We establish
theoretical guarantees for density estimation with fully Bayesian Markoy models (Section 3)
AR models used for compression, like other AR models, can be embedded into BEAR models for improved statistical performance and to measure misspecification.
BEAR models are closely linked to non-generative genome analysis methods. Assembly algorithms and variant callers often analyze paths in the de Bruijn graph of sequence dataset; in the limit h + oc, samples from the posterior predictive distribution of the BEAR model, conditional on L, correspond to paths through the L-mer de Bruijn graph of the data I 33 Comparisons between genomes and other sequences are often made on the basis of kmer counts; our two-sample test provides generative perspective on this idea R 16
BEAR models are also connected to ideas in natural language processing, where kmers are referred to ngrams. Under the vanilla. BEAR model, the mean of the posterior predictive distribution conditional on L corresponds to an ngram additive smoothing model M. Comparisons between datasets using their ngram counts are also common in model evaluation metrics such as the BLEU score
Density estimation
The density estimation problem is that of estimating p* given data (X,) drawn i.i.d. from p". Density estimation is particularly crucial for biological sequence analysis due to its connections to fitness estimation 31 N5 Slate-of-the-art mutation effect prediction methods and clinical variant interpretation methods rely on density est istattes of evolutionary sequence data 20 M. Density estimation with generative models is particularly useful for protein design, as samples from accurate density estimates likely to be functional and

---

Despite all these applications, existing density estimation methods for biological sequences lack theoretical guarantees on their accuracy and are often limited in their scale, being restricted to relatively short sequences [ Here, we show that the posterior distribution of the BEAR model is consistent and will concentrate on p" as N + co, regardless of what p" actually is, so long as p" generates finite length sequences almost surely (a.s.).
can be synthesized in the laboratory 83 66
We first study the expressiveness of BEAR models. Let Mr be the set of Markov models Py with transition probabilities a and lag L that generate finite length strings a.s.. Note that M1 C M2 .... Define the union M = Up=1 We can compare M to the set of distributions over finite strings S. of which p" is a member. In Section D we prove that, Summary of Propositions 14 Not all possible distributions over S are in M. However, M is dense on the space of probability distributions over S with the total variation metric. The implication of this result is that although BEAR models cannot exactly match arbitrary lata-generating distributions, they can approximate p" arbitrarily well as increases. This makes asymptotic consistency possible.
We now show that the posterior of the BEAR will in fact asymptotically concentrate on the true p", i.e. it is consistent. For tractability, we assume in this section that the prior is fixed (we do not use empirical Bayes). The result relies on the tools for understanding convergence rates of posteriors developed in Ghosal et al. [2 The most important assumption is that p" is subexponential, meaning that for some > O, Ep- exp(t)X() < co where |XX] is the sequence length. Let [1(II(X,)m=z)) denote the posterior over sequence distributions. Let B(p". denote a ball of radius 5 centered at p", using the Hellinger distance.
Summary of Thcorem 35 Given M > O large enough and E (0, 1) small enough, we
have [1(" MN 1-) (X,) =)===)>>>> +
in probability.
A proof is in Section H and simulations in Section 1.2 This result states that the posterior
distribution of the model converges to a delta function at the true distribution p' regardless
of what p" is. It also provides rate of convergence: in a parametric model, the uncertainty would shrink as N-1, but here the rate is slower, - price paid for the nonparametric model's expressivity 23 2g 35. The proof includes a variety of new theoretical constructions and algorithms that are used to approximate subexponential sequence distributions.
4 Robust parameter estimation
To derive a biological understanding of mutational processes, evolutionary history, funational constraints, etc. from sequence data, researchers must estimate model parameters (not just density). However, parameter estimates cannot in general be trusted when models are misspecified [3 To reach robust scientific conclusions, therefore, parameter estimates should idcally come with warning about whether or not the model is misspecified and some measurement of the degree of misspecification. Here, we study in BEAR models the asymptotic behavior of empirical Bayes estimates of the AR parameter 8, as well as the hyperparameter h, showing that h diagnoses misspecification in the embedded AR model.
Our analysis builds off the study of empirical Bayes consistency in Petrone et al. bd which showed that empirical Bayes will, in general, maximize the prior probability of the true data-generating parameter value. Extending this theory to BEAR models is nontrivial, since in BEAR models the standard Laplace approximation to the marginal likelihood can fail. For theoretical tractability, as in many analyses of similar models, we fix L at some ar arbitrary and large value B0]. Define p (L) = argmin,, E.EM, KL(p" (Ip.)) as the closest model in M, to p", and define e" such that P.- = (L)) (note (=(() - p" as - oo). We say that the AR model is misspecified "at resolution L" if f cannot approximate p=(L), Le. if there does not exist some sequence of parameter values ON such that Pς(āw) + +p)(2) as N + oo; otherwise, the AR model is well-specified at resolution L. Now we can study empirical Bayes estimates of h and 0, denoted fNn and ON
Let ( and (0 ( be sequences mazimizing the BEAR marginal likelihood pl ((((x)i w L,, h 0) for each N. If the model is well-specified at resolution I., then by N1/1-e - for every 6 > O and Ps(0x) -> p (L(L) in distribution, with both sequences converging in probability. On the other hand, if the model is misspecified at
Summary of Propositions 15 20
5

---

resolution L, then hy is eventually bounded below by some positive (non-zero) number a.s.. Proofs are in Section and simulations in Section The implication of this result that when the AR model is well-speciffed, hy converges to zero (at rate that is a power of the dataset size) and 0 converges to the parameter value 0* at which the AR model matches the data (Corollary T6.. On the other hand, when the AR model is misspecified, hy does not converge to zero; heuristically, we find instead that by is approximately proportional to divergence between p²(4) and the AR model,
Σ (garoz( (pr)
(L(f((ON)ll(t) + (og((N)
Σ
fas(eN))
(4)
besuppe (r")).
where accz = [k (*k>> 0} is the set of kmers with non-zero probability and (upp;("")) 3|2 = {6 p" (k), > 0) > 0} is the set of transitions from k with non-zero probability. In summary: when fitting a BEAR model by empirical Bayes, you get, along with parameter estimate ON, value hy which tells you the amount (from zero to infinity) of misspecification in the AR model. If hy is close to zero, you can trust the estimate ON-
Hypothesis testing
Goodness-of-fit test A major outstanding challenge in biological sequence analysis is to build models based on natural sequence data that are accurate enough to generate novel functional sequences [5. A crucial component of the problem is model evaluation: while relative model performance may be compared on the basis of likelihood, absolute performance - whether or not the model in fact provides an accurate description of the data - is usually addressed solely on the basis of limited numbers of summary statisties, such as average amino acid hydrophobicity or sequence length [R Bitl Given a datasel (Xa) p" i.i.d., goodness-of-fi test asks whether not the data distribution p" matches a model distribution P. It takes into account all possible distributions p" including those that differ from ò in a manner that cannot be captured by finitely many summary statistics. We propose a goodness-of-fit test that compares the mull hypothesis Ho : p" P to the alternative H, :p" + P using the Bayes factor BF = >>((X,) h, 0)/p())(p))a (X1m), where IL, h, ()T(() L) is the marginal likelihood under the BEAR model. Note that practically, the sum over L straightforward to approximate by truncation, and that the test can be computed in time linear in the amount of data.
((xx)))__ h, 0) = Ez p((x))
We now prove the consistency of the test. As in comparable theoretical analyses of tests based on Polya trees, for theoretical tractability we truncate the prior, setting x = 0 for L larger than some arbitrary L but (L) > 0 for L S Î 30. We treat O and h > 0 as fixed. Summary of Proposition 21 Iʃ p is al least as close to p' as p=(L) is, as measured by kl.(""ll),, then BF - O in probability as N + co. On the other hand, if p<(L) is closer than D, then BF + oo in probubility. A proof is in Section G.T and simulations in Section L3
An important practical limitation nonparametrie hypothesis testing is low power: since so many alternative distributions must be considered, the null hypothesis can rarely be rejected. However, Proposition 21 holds for the Bayes factor (F,,, 0) p((X)) ) L, h, ())(((()))))))))))))))))(())() (=1) with any choice of L, h > 0, and 8. Thus in practice to increase power we can maximize the value of (F h, 0) as function of L, h, and/or 0 (note that this approach is heuristic, since we have not proven the consistency of the maximized Bayes factor). Berger and Guglielmi E provide extensive methodological guidance analogous tests constructed with Polya trees. Based on their recommendations we suggest first choosing 0 such that Py(0) as close possible to p, then plotting the Bayes factor as function of l and/or L to identify the maximum value and confirm that any conclusion is robust to changes in h and/or orr L.
using
Another challenge in nonparametric hypothesis testing is that it can be difficult to understand how exactly a test reached its conclusion. To identify which sequences provided the most evidence for or against the null hypothesis, we suggest examining the BEAR Bayes factor for each individual sequence conditional on the rest of the dataset, in analogy to the witness function used in kernel-based tests [3 70

---

Table 1: Heldout perplexity. Whole genome sequencing data: YSD1: A Salmonella phage. A. the: Arabidopsis thaliana, a plant (datasets represent different individuals). Single cell RNA sequencing data: PBMC: peripheral blood mononuclear cells, taken from a healthy donor. HL: Hodgkin's lymphoma tumor cells. GBM: glioblastoma tumor cells. Metagenomic sequencing data: HC: non-CD and non-UC controls. CD: Crohn's disease. UC: ulcerative colitis. Full assembled genomes: Bact.: Bacteria, Models Van Vanilla (Jeffreys prior). Lin.:
Linear. CNN: convolutional neural network. Ref.: reference genome/ wwwaeeertatnnntnneeeeeeeee model.
Dataset AR Lin.
AR CNN
AR Ref.
BEAR Van. BEAR Lin. BEAR CNN BEAR Ref.
YSD1
3.953
3.873
1.266
1.165
1.144
1.144
1.145
A. th. 1
3.956
3.947
2.686
1.567
1.432
1.432
1,411
A. th. 2
3.953
3.949
1.982
1.650
1.463
1.462
1.441
A. th. 3
3.998
3.952
2.340
1.834
1.728
1.727
1.733
PBMC
3.991
3.974
2.097
1.402
1.372
1.372
1.374
HI
3.959
3.930
2.141
1.409
1.378
1.378
1.379
GBM
4.137
4.137
2.366
1.442
1.406
1.406
1.406
HC
3.966
3.946
1.652
1.465
1.464
CD
3.992
3.985
1.760
1.524
1.524
UC
3.989
3.986
1.644
1.481
1.481
Bact.
3.831
3.794
3.774
3.774
3.774
Two-sample test two-sample test asks whether or not datasots (X,) and (X,) y are drawn from the same distribution. Efforts to compare different sequence datasets are widespread in biology: for instance, researchers often wish to determine whether two microbiome samples, taken under different conditions at different timepoints, are the same up to sampling noise [d Two-sample tests can also be used to evaluate generative sequence models that lack tractable likelihoods (for which the goodness-offff of-fft test proposed above does not apply) such as energy- anergy-based models or implicit models like GANs and biophysical simulators p7 12 49 Assume (Xn) ~ Pı and (X,) ~ i.i.d.. Our BEAR test compares the null hypothesis Ho P1 = P2 the alternative Hi # p2 using the Bayes factor BF = p((((((,))) (= ()x((X;;) "" 0) (XX)) (X, (X,)==/hh, 0). As in the goodness-off tfaness-offitt case, the test can be computed approximately in time linear in the amount of data, and the same advice on increasing power and identifying important sequences holds here too.
We next prove consistency, again truncating the prior at L and fixing h and A.
Summary of Proposition 22 If pl Ps (f.) then BF + 0 as N + oo in probability. (f. Otherwise, if (f) +
(i.
then BF - oc in probability. proof is in Section G.2 and simulations in Section
Results
Predicting sequences We sought to evaluate BEAR models as compared to AR models on the task of predicting real nueleotide (nt) sequences. We considered eleven datasets of four different types: whole genome sequencing read data, single cell RNA sequencing read data (including from patient tumors), metagenomic sequencing read data (including from patient fecal samples) and full bacterial genomes from across the tree of life (Section K Datasets ranged in total size from 107 10a0 and in individual sequence length from 102 10 nt (Table S1). 25% of data was randomly held out for testing, in the form of entire sequences (reads, genomes, etc., see Table 52): our goal was to evaluate BEAR models as density estimators, so we did not use masking (a common holdout strategy in natural language processing). We considered linear AR model and deep convolutional neural network (CNN) AR model with > 10x more paramelers, both of which are common models used across a range of applications; we also designed biologically-strudarres AR model which makes predictions based reference genome and Jukes- Cantor mutation model

---

(Section L.1
56 37. We then embedded each AR model to create a corresponding BEAR model. The BEAR models improve over the AR models in nucleotide prediction according to both perplexity (Table I and accuracy (Table S3 in all datasets, even when the model lag I is held fixed for comparison (Section
In 10 out of 11 datasets, BEAR models increase nucleotide prediction accuracy from near chance values of 30 35% (in the case of the linear and CNN models) to 78 95%, bringing genome-scale models into the realm of potential practical use (Table S3 The training time for BEAR models is essentially identical to that of AR models, aside from the time required to build the transition count matrix, which need only be done once before training all models (Fig. (SI).. Remarkably, the optimal lag L chosen by empirical Bayes is often quite short, less than 20 nt (Table S4. The improvements offered by BEAR models that use an embedded AR model over the vanilla BEAR model are modest for datasets of this size; however, sequencing experiments are often designed to collect enough data for downstream analyses. We found in an example that, if sequencing coverage was 3x instead of 100x, the improvement in prediction accuracy would have been greater than 10 percentage points instead of 0.1 (Section L.4 Fig. S14..
Measuring misspecification When conventional deep neural network methods fail to provide strong predictive performance, popular wisdom often ascribes the failure to too much model flexibility or not enough training data, especially in scientific applications. Ex- amining the h misspecification diagnostic in the BEAR models described above, we see that this is not the case here (Table 2 The large values of suggest that where the CNN fails it is not because of too much flexibility but rather too little: the model is not flexible enough to encompass the true data distribution, so it suffers from misspecification. Meanwhile, the reference- based model has only two learned parameters, but is less misspecified than the CNN in all but one dataset. This too runs counter to popular wisdom in machine learning, which often assumes that when principled, low-flexibility sci- entific models outperform deep neural networks it is thanks to their low variance in the small data regime.
Table 2: Diagnostic h. Abbrevia- tions in Table II]
Dataset Lin. CNN
Ref.
YSD1 5.528 5.461 4.183
A. th. 1 2.765
2.756
2.990
A. th. 2.643
2.633
2.326
A. th. 3 3.969
3.964
1.598
PBMC
4.167
4.145
3.762
HL
4.050
4.038
3.581
GBM
4.172
4.154
3.238
HC
4.668
4.651
CD
3.096
3.094
UC
3.843
3.835
Bact.
0.010
0.003
Generating samples BEAR models are generative and can be used to sample new sequences. We sampled extrapolations from the end of a read sequence recorded in a plant (A. thaliana) whole genome sequencing experiment, and compared to alternative non-probabilistic extrapolation method that is widely used in biology, local assembly (Fig. A Section M.. In this example the assembly algorithm SPAdes returns four possible assemblies. relatively large mumber compared to other reads in the dataset (Fig. 3 stars) F Samples from the BEAR model include these four possibilities, but also many more, some with higher probability. The distribution over possible mucleotide choices under the BEAR model is much wider than the number of assemblies would suggest: it has perplexity of 1.4 per position (on average across samples) at the beginning of the extrapolation, and a perplexity of 2.7 at 50 nucleotides (Fig- BB). These observations suggest that SPAdes, which does not provide measurement of uncertainly, may not be capturing the full range of possible sequences.
Visualizing data Methods for learning local representations features of biological se- quences can be powerful tools for visualization and semisupervised learning M. One approach to extracting such representations is to learn a generative model qX1...... XL-1) 1) of kmers, for instance using a variational autoencoder. While such models are not autorogressive, the small size of the DNA alphabet makes it tractable to estimate the conditional ((x2+11X11) by Bayes' rule, and this conditional can then be embedded into a BEAR model. We applied this strategy to probabilistie PCA. We visualized in low dimensions the inferred latent representation for a model trained on a single cell RNA sequencing dataset (HL), and were able to assign annotations to clusters, including those containing unmapped reads (Fig. BCD; Section N). The BEAR model however raises the warning that the model is misspecified (h 4.836), suggesting there may be richer latent structure yet to discover.
8

---

log P
A. Generated samples
B. Sample diversity
C. Prediction comparison
D. Latent representations
-1
3.1
Reference model log
-88
Unilgred
433
-10
-33
perpleaty
TSNE 2
-120
O Miscchnndia
10
1000
O APAO
also SPAdes local assemblies
position
BEAR model log
TSNE
E. Goodness-offfit tests
F. Two-sample tests
G. Evidence in two-sample tests
A. thaliana
A. thaliane
2,82
103
(mismatch/mate m
log Bayes factor
2 a
log Bayes factor (mismatch/match)
Number of reads
bactonie' DETALE
2.24
-52
52 LCD 152
10%10
Lag
log likelihoad ratio (real data/simulated)
Figure 3: Generation, visualization and testing. (A) Sample extrapolations, colored to denote distinct paths through the L-mer de Bruijn graph. (B) Distribution of the perplexity of the next Markov transition under the BEAR model, for each position of the sampled extrapolations, with the per position average shown in black (Section Mp. (C) Log probability
of each read in the HL dataset under the BEAR model and model built from the reference transcriptome. Reads are colored by whether or not they map to the reference. (D) Latent reprosentations of the reads highlighted in C, visualized using tSNE, with clusters annotated as likely coming from mitochondria, the sequencing adapter, or transcripts of the gene JUND (Section No. (E) ofo-nt test Bayes factor as function of hy pyrraarraeteer h. (F) Two-sample test Bayes factor as a function of lag L. Black line compares simulated data to simulated data; dashed lines compare subsampled real data to subsampled real data; solid lines compare real data to simulated data. (G) Log probability of each read under the real data BEAR model minus the log probability under the simulated data BEAR model (Section O
Testing hypotheses The question of when and how microbiomes change is widespread, but has in the past relied on summary statistics of sequencing datasets μل]. Schreiber et al. [o] studied changes in patient urine microbiomes before and after kidney transplant, and performed both unbiased melagenomic sequencing and diagnostic quantitative ive polymerase chain reaction (qPCR) for a specific virus associated with complications (JC polyomavirus).
They found evidence of tonor-to-recipient viral transmission in 5 cases out of 14. We applied the BEAR two-sample test to patients' metagenomic sequencing data before and after transplantation, using the vanilla Jeffreys prior and integrating over lags, in order to detect changes; the test rejects the null hypothesis in all 5 cases where there was transmission, and accepts the null hypothesis in all but one of the remaining 9 cases (Table St Section O.I. These results show, in small example, that the two-sample test has sufheient power to delect microbiome changes in real data, and can be consistent with more specific lests.
We next applied BEAR hypothesis tests evaluate generative models. We evaluated the reference-based AR model described above using the BEAR goodness-of-fit test. The test identifies considerable evidence (log Bayes factor > 10a) for misspecification in each A. thaliana whole genome sequencing dataset, and this conclusion is robust to a wide range of values (Fig. c Section 6.2). Next, we evaluated a detailed simulation model (ART) that is intended to generate likely reads of a given reference genome 32 The model lacks tractable likelihoods, so we use the BEAR two-sample test. When integrating over all lags, the test accepts the null hypothesis, suggesting that the simulation model is accurate; if we examine the test results for individual lags L to increase power, however, we can see some evidence of differences (Fig. 3 Section 0.2 Note that as L increases, there is a tradeoff: tests with larger lag can detect more subtle differences between the two distributions, but have less statistical power since they must consider larger set of possible distributions. Thus the
9

---

Bayes factor first increases and then decreases with lag, reaching peak at intermediate values where there is the most evidence of difference. To understand in detail the source of the detected differences between the data and the simulation model, we examined the conditional Bayes factor for individual reads, discovering clusters of reads that are poorly explained by the simulation model (Fig. G One group mapped to chloroplasts, with its own genome that is variable in copy number; reads mapping to centromeres, an area of the plant genome for which the reference genome is considered unreliable, were also poorly explained by the simulation model. In one dataset we found a cluster of outliers that did not map to A. thaliana at all, and instead mapped to common soil bacteria, Bacillus cereus, presumably a contaminant in the experiment (Fig. 3G, left). These results illustrate how BEAR hypothesis tests can be used not only for testing but also for detailed model criticism.
organelle
7
Discussion
In this article we proposed the nonparametrie BEAR model, studied its theoretical properties, and developed algorithms and implementations for terabyte-scale inference. BEAR models substantially outperform standard AR models on a variety of datasets, and come with extensive theoretical guarantees, including for density estimation, misspecification detection, and hypothesis testing. BEAR models are closely connected to non-probabbllitic genome analysis methods, such as de Bruijn graph assembly, but provide an alternative that is uncertainty-aware. Note, however, that BEAR models do not explicitly account for paired- end read information, or other sources of long-distance information; this is an important area for future work. While there has been little previous empirical or theoretical work in the machine learning literature on generative models of full genomic, transcriptomic or metagenomic sequences, we hope BEAR models provide a useful starting point.
Acknowledgments and Disclosure of Funding
We thank Jean Disset for a small scale version of the kmer counting code, Rob Patro for crucial advice on large scale kmer counting, and Winnie Wang for illustrations used in the theory section of the supplementary material. We also thank Chris Sander, Elizabeth Wood, Tessa Green and all the members of the Marks Lab for discussion and ideas. E.N.W. is supported by the Fannie and John Hertz Foundation. D.S.M. is supported by the Chan Zuckerberg Initiative.
References
[1] 1001 Genomes Consortium. 1,135 genomes reveal the global pattern of polymorphism in Arabidopsis thaliana. Cell, (26(2))881 491, July 2016.
[2] M. Abadi, P. Barham, J. Chen, Z. Chen, A. Davis, J. Dean, M. Devin, S. Ghemawat, G. Irving, M. Isard, and Others. Tensorflow: A system for large-scale machine learning. In 12th USENIX symposium operating systems design and implementation (OSDI
16), pages 265-26250 283. usenix.org, 2016.
[3] B. Alsop and J. Raymond. Resolving Prokaryotic Taxonomy without rRNA: Longer Oligonucleotide Word Lengths Improve Genome and Metagenome Taxonomic Classifica- tion. PLoS ONE, 8(7), 2013.
[4] J. L. Ba, J. R. Kiros, and G. E. Hinton. Layer normalization. July 2016.
[5] A. Bankevich, S. Nurk, D. Antipov, A. A. Gurevich, M. Dvorkin, A. S. Kulikov, V. M. Lesin, S. I. Nikolenko, S. Pham, A. D. Prjibelski, A. V. Pyshkin, A. V. Sirotkin, N. Vyahhi, G. Tesler, M. A. Alekseyev, and P. A. Pevzner. SPAdes: new genome assem- bly algorithm and its applications to single-cell sequencing. Journal of Computational Biology, 4577 May 2012.
[6] J. O. Berger and A. Guglielmi. Bayesian and conditional frequentist testing of parametric model versus nonparametric alternatives. Journal of the American Statistical Association, 96(453):174-188) 2001.
10

---

[7] S. Biswas, G. Khimulya, E. C. Alley, K. M. Esvelt, and G. M. Church. Low-N protein engineering with data-efficient deep learning. Nature Methods, 18(4):389 396, 2021.
[8] G. M. Boratyn, C. Camacho, P. Cooper, G. Coulouris, A. Fong, N. Ma, T. L. Madden, W. T. Matten, S. D. McGinnis, Y. Merezhuk, Y. Raytselis, E. W. Sayers, T. Tao, J. Ye, and I. Zaretskaya. BLAST: a more elficient report with usability improvements. Nucleic Acids Research, 41(Web Server issue): W29 33, July 2013.
99] S. F. Chen and J. Goodman. An empirical study of smoothing techniques for language modeling. Computer Speech and Language, 13(4):359 394, Oct. 1999.
[10] P. J. A. Cock, T. Antao, J. Chang, B. Chapman, C. Cox, A. Dalke, I. Friedberg, T. Hamelryck, F. Kauff, B. Wilezynski, and M. J. L. de Hoon. Biopython: freely available python tools for computational molecular biology and bioinformatics. Bioinformatics, 25(11):1422 1423, June 2009.
[11] P. E. C. Compeau, P. A. Pevzner, and G. Tesler. How to apply de Bruijn graphs to genome assembly. Nature Biotechnology, 29(11):987799) 991, Nov. 2011.
[12] I. Cortés-Ciriano, J. J.-K. Lee, R. Xi, D. Jain, Y. L. Jung, Yang, D. Gordenin, L. J. Klimezak, C.-Z. Zhang, D. S. Pellman, PCAWG Structural Variation Working Group, P. J. Park, and PCAWG Consortium. Comprehensive analysis of chromothripsis in 2,658 human cancers using rrhole-geone sequencing. Nature Genetics, 52(3):3311 341, Mar. 2020.
[13] S. Darmanis, S. A. Sloan, D. Croote, M. Mignardi, S. Chernikova, P. Samghababi, Y. Zhang, N. Neff, M. Kowarsky, C. Caneda, G. Li, S. D. Chang, I. D. Connolly, Y. Li, B. A. Barres, M. H. Gephart, and S. R. Quake. Single-Cell RNA-Seq analysis of infiltrating neoplastic cells at the migrating front of human glioblastoma. Cell Reports, 21 (1):11339911444 Oct. 2017.
[14] A. P. Dawid. Posterior model probabilities. In P. S. Bandyopadhyay and M. R. Forster, editors, Philosophy of Statisties, volume 7, pages 607-630. North-Holland, Amsterdam, Jan. 2011.
[15] J. V. Dillon, I. Langmore, D. Tran, E. Brevdo, S. Vasudevan, D. Moore, B. Patton, A. Alemi, M. Hoffman, and R. A. Saurous. TensorFlow distributions. Nov. 2017.
[16] V. B. Dubinkina, D. S. Ischenko, V. I. Ulyantsev, A. V. Tyakht, and D. G. Alexeev. Assessment of k-mer spectrum applicability for metagenomic dissimilarity analysis. BMC Bioinformatics, 17:38, Jan. 2016.
[17] M. Duc Cao, T. Dix, L. Allison, and C. Mears. A simple statistical algorithm for biological sequence compression. In 2007 Data Compression Conference (DCC'07), pages 43 52, Mar. 2007.
[18] R. A. Dunstan, D. Pickard, S. Dougan, D. Goulding, C. Cormie, J. Hardy, F. Li, R. Grin- ter, K. Harcourt, L. Yu, J. Song, F. Schreiber, J. Choudhary, S. Clare, F. Coulibaly, R. A. Strugnell, G. Dougan, and T. Lithgow. The flagellotropic bacteriophage YSD1 targets salmonella typhi with a chi-like protein tail fibre. Molecular Microbiology, 112 16):1831-1866, Dec. 2019.
[19] R. Durbin, S. Eddy, A. Krogh, and A. Mitchison. Biological Sequence Analysis. 1998.
[20] J. Frazer, P. Notin, M. Dias, A. Gomez, K. Brock, Y. Gal, and D. Marks. Large-scale clinical interpretation of genetic variants using evolutionary data and deep learning. Dec. 2020.
[21] C. Galiez, M. Siebert, F. Enault, J. Vincent, and J. Söding. WIsH: who is the host? predicting prokaryotic hosts from metagenomic phage contigs. Bioinformatics, 33(19): 3113 3114, Oct 2017.
[22] S. Geman and C.-R. Hwang. Nonparametric maximum likelihood estimatiin by the method of sieves. The Annals of Statistics, 10(2):401-411 114, 1982.
11

---

[23] S. Ghosal, J. K. Ghosh, and A. W. van der Vaart. Convergence rates of posterior distributions. The Annals of Statistics, 28(2):500-555) 2000.
[24]] J. Ghosh and R. Ramamoorthi. Bayesian Nonparamotries. 2003.
[25] P. Gopalan, W. Hao, D. M. Blei, and J. D. Storey. Scaling probabilistic models of genetic variation to millions of humans. Nature Genetics, 48(12):1587-155) 1590,, Dec. 2016.
[26] R. M. Gray, Entropy and Information Theory. Springer Science & Business Media, Jan. 2011.
[27] A. Gretton, K. M. Borgwardt, M. J. Rasch, B. Schölkopf, and A. Smola. A kernel two-sample test. Journal of Machine Learning Research, 13:723-773, 2012.
[28] C. R. Harris, K. J. Millman, S. J. van der Walt, R. Gommers, P. Virtanen, D. Cournapeau,
E. Wieser, J. Taylor, S. Berg, N. J. Smith, R. Kern, M. Picus, S. Hoyer, M. H. van Kerkwijk, M. Brett, A. Haldane, J. F. del Río, M. Wiebe, P. Peterson, P. Gérard- Marchant, K. Sheppard, T. Reddy, W. Wockesscr, H. Abbasi, C. Gohlke, and T. E. Oliphant. Array programming with NumPy. Nature, 585 (85(7825):3577338
Sept. 2020.
[29] T. Henighan, J. Kaplan, M. Katz, M. Chen, C. Hesse, J. Jackson, H. Jun, T. B. Brown, P. Dhariwal, S. Gray, C. Hallacy, B. Mann, A. Radford, A. Ramesh, N. Ryder, D. M. Ziegler, J. Schulman, D. Amodei, and S. McCandlish. Scaling laws for autoregressive gencrative modcling. Oct. 2020.
[30] C. C. Holmes, F. Caron, J. E. Griffin, and D. A. Stephens. Two-sample bayesian nonparametric hypothesis testing. Bayesian Anal., (((2):2773200) June 2015.
[31] T. A. Hopf, J. B. Ingraham, F. J. Poelwijk, C. P. I. Schärfe, M. Springer, C. Sander, and D. S. Marks. Mutation effects predicted from sequence co-variation. Nature Biotechnology, 35(2):128 135, Feb. 2017.
[32] W. Huang, L. Li, J. R. Myers, and G. T. Marth. ART: next-gencration sequencing read simulator. Bioinformatics, 28(4):593 594, Feb. 2012.
[33] 7. Iqbal, M. Caccamo, I. Turner, P. Flicek, and G. McVean. De novo assembly and genotyping of variants using colored de Bruijn graphs. Nature Genetics, 44(2):226-232, Jan. 2012.
(34] P. E. Jacob, L. M. Murray, C. C. Holmes, and C. P. Robert. Better together? statistical learning in models made of modules. Aug. 2017.
[35] J. Kaplan, S. McCandlish, T. Henighan, T. B. Brown, B. Chess, R. Child, S. Gray, Radford, Wu, and D. Amodei. Scaling laws for neural language models. Jan. 2020.
(36] D. Kim, J. M. Paggi, C. Park, C. Bennett, and S. L. Salzberg. Graph-based genome aligmment and genotyping with HISAT2 and HISAT-genotype. Nature Biotcchnology, 37(8):907 915, Aug. 2019.
[37] D. P. Kingma and J. Ba. Adam: A method for stochastic optimization. In International Conference on Learning Representations, 2015.
[38] D. P. Kingma and M. Welling, Aulo-encoding variational bayes. In International Conference on Learning Representations. 2014.
[39] M. Kokot, M. Dlugosz, and Deorowicz. KMC 3: counting and manipulating k-mer statistics. Bioinformatics, 33(1))22759 2761, Sept. 2017.
[40] W. Kool, H. van Hoof, and M. Welling. Stochastic beams and where to find them: The Gumbel-Top-2 trick for sampling sequences without replacement. In International Conference on Machine Learning, pages 3499 3508. PMLR, 2019.
[41] A. Kucukelbir, D. Tran, R. Ranganath, A. Gelman, and D. M. Blei. Automatic differentiation variational inference. Journal of Machine Learning Research, 18(14): 1-45, Jan. 2017.
12

---

[42] Y. Li, K. Swersky, and R. Zemel. Generative moment matching networks. In International Conference on Machine Learning, pages
PMLR, 2015.
[43] J. R. Lloyd and Z. Ghahramani. Statistical model criticism using kernel two sample tests. In Advances in Neural Information Processing Systems, pages 829- 837, 2015.
[44] J. Lloyd-Price, A. Mahurkar, G. Rahnavard, J. Crabtree, J. Orvis, A. B. Hall, A. Brady, H. H. Creasy, C. McCracken, M. G. Giglio, D. McDonald, E. A. Franzosa, R. Knight, O. White, and C. Huttenhower. Strains, functions and dynamics in the expanded human microbiome project. Nature, 550(7674):611666) 66, Oct. 2017.
[45] A. Madani, B. Krause, E. R. Greene, S. Subramanian, B. P. Mohr, J. M. Holton, J. L. Olmos, C. Xiong, 7. Z. Sun, R. Socher, J. S. Fraser, and N. Naik. Deep neural language modeling enables functional protein generation across families. 2021.
[46] G. Marçais and C. Kingsford. A fast, lock-free approach for efficient, parallel counting of occurrences of k-mers. Bioinformatics, 27(())7644-77 770, Mar. 2011.
[47] J. W. Miller. Asymptotic normality, concentration, and coverage of generalized posteriors. Journal of Machine Learning Research, 22(168):1-55 53, 2021.
[48] J. W. Miller and D. B. Dunson. Robust bayesian inference via coarsening. Journal of the American Statistical Association, 114(527):1113-111 2019.
[49] S. Mohamed and B. Lakshminarayanan. Learning in implicit generative models. Oct. 2016.
[50] R. E. Mukamel, R. E. Handsaker, M. A. Sherman, A. R. Barton, Y. Zheng, S. A. McCarroll, and P.-R. Loh. Protein-coding repeat polymorphisms strongly shape diverse human phenotypes. Science, (7362)) 14999 1499-1505, Sept. 2021.
[51] Narlikar, N. Mehla, S. Galande, and M. Arjunwadkar. One size does not fil all: on how Markov model order dictates performance of genomic sequence analyses. Nucleic Acids Research, (33:1441 1424, Feb. 2013.
[52] N. A. O'Leary, M. W. Wright, J. R. Brister, S. Ciufo, D. Haddad, R. McVeigh, B. Rajput,
B. Robbertse, B. Smith-White, D. Ako-Adjei, A. Astashyn, A. Badreldin, Y. Bao, O. Blinkova, V. Brover, V. Cheivernin, J. Choi, E. Cox, O. Ermolaeva, C. M. Farrell, Goldfarb, T. Gupta, D. Haft, F. Hatcher, W. Hlavina, S. Joardar, V. K. Kodali, W. Li, D. Maglott, P. Masterson, K. M. McGarvey, M. R. Murphy, K. O'Neill, S. Pujar, S. H. Rangwala, D. Rausch, L. D. Riddick, C. Schoch, A. Shkeda, F. Thibaud-Nissen, 1. Tolstoy, R. E. Tully, A. R. Vatsan, C. Wallin, D. Webb, W. Wu, M. J. Landrum, A. Kimchi, T. Tatusova, M. DiCuccio, P. Kitts, T. D. Murphy, and K. D. Pruitt. Reference sequence (RefSeq) database at NCBI: current status, taxonomic expansion, and functional annotation. Nucleic Acids Research, 44(D1))D733 45, Jan. 2016.
S. Storz, H. Sun,
[53] K. Papineni, S. Roukos, T. Ward, and W.-J Zhu. BLEU: a method for automatic evalu- ation of machine translation. In Anmual Meeting of the Association for Computational Linguistics, pages 311-318, 2002.
[54] S. Petronee J. Rousseau, and C. Scricciolo. Bayes and empirical Bayes: do they merge? Biometrika, 101(2):285 302, 2014.
[55] A. J. Pinho, P. J. S. G. Ferreira, A. J. R. Neves, and C. A. C. Bastos. On the representability of complete genomes by multiple competing finite-context. (Markov) models. PLoS One, 6(6):e2158) 5888, June 2011.
[56] D. Pratas, M. Hosseini, and A. J. Pinho. Substitutional tolerant markov models for relative compression of DNA sequences. In International Conference on Practical Applications of Computational Biology & Bioinformatics (PACBB), pages 265-22 (65-272,
June 2017.
13

---

[57] J. K. Pritchard, M. Stephens, and P. Donnelly. Inference of population structure using multilocus genotype data. Geneties, 153(2):945 959, June 2000.
[58] J. Ren, X. Bai, Y. Y. Lu, K. Tang, Y. Wang, G. Reinert, and F. Sun. Alignment-free sequence analysis and applications. Annual Review of Biomedical Data Science, 1: 114, July 2018.
[59] D. J. Rezende, S. Mohamed, and D. Wierstra. Stochastic backpropagation and ap- proximate inference in doep generative models. In International Conference on Machine Learning, 2014.
[60] A. J. Riessclman, J. B. Ingraham, and D. S. Marks. Doep generative models of genetic variation capture the effects of mutations. Nat. Methods, 15(10):816 822, Oct. 2018.
[61] J. Rousscau. On the frequentist properties of Bayesian nonparametric methods. Annual Review of Statistics and Its Application, 3:211-231, 2016.
[62] J. Rousseau and K. Mengersen. Asymptotic behaviour of the posterior distribution in overfitted mixture models. Journal of the Royal Statistical Society: Series B, 73(5): 189-710, Nov. 2011.
[63] W. P. Russ, M. Figliuzzi, C. Stocker, P. Barrat-Charlaix, M. Socolich, P. Kast, D. Hilvert, R. Monasson, S. Cocco, M. Weigt, and R. Ranganathan. An evolution-based model for designing chorismate mutase enzymes. Science, 369:440 445, 2020.
(64] P. W. Schreiber, V. Kufner, K. Hübel, S. Schmutz, O. Zagordi, A. Kaur, C. Bayard, M. Greiner, A. Zbinden, R. Capaul, J. Böni, H. H. Hirsch, T. F. Mueller, N. J. Mueller, A. Trkola, and M. Huber. Metagenomic virome sequencing in living donor and recipient kidney transplant pairs revealed JC polyomavirus transmission. Clinical Infectious Diseases, (6(6):987 994, Aug, 2019.
[65] G. Sella and A. E. Hirsh. The application of statistical physics to evolutionary biology. Proceedings of the National Academy of Sciences, 102(27):9541 9546, July 2005.
[66] J.-E. Shin, A. J. Riesselman, A. W. Kollasch, C. McMahon, E. Simon, C. Sander, A. Manglik, A. C. Kruse, and D. S. Marks. Protein design and variant prediction using autoregressive generative models. Nature Communications, 12(1):2403, Apr. 2021.
[67] M. Silva, D. Pratas, and A. J. Pinho. Efficient DNA sequence compression with neural networks. Gigascience, ((1), Nov. 2020.
[68] J. T. Simpson, K. Wong, S. D. Jackman, E. Schein, S. J. M. Jones, and Birol. ABySS: a parallel assembler for short read sequence data. Genome Research, 19(6): 1117-1123 June 2009.
[69] P. J. Stephens, C. D. Greenman, B. Fu, F. Yang, G. R. Bignell, L. J. Mudie, E. D. Pleasance, K. W. Lau, D. Beare, L. A. Stebbings, S. McLaren, M.-L. Lin, D. J. McBride, I. Varela, S. Nik- Nik-Zainal, C. Leroy, M. Jia, A. Menzies, A. P. Butler, J. W. Teague, M. A. Quail, J. Burton, H. Swerdlow, N. P. Carter, L. A. Morsberger, C. Tacobuzio-Donahus G. A. Follows, A. R. Green, A. M. Flanagan, M. R. Stratton, P. A. Futreal, and P. J. Campbell. Massive genomic rearrangement acquired in single catastrophic event during cancer development. Cell, 144(1):27-40, Jan. 2011.
(70] D. J. Sutherland, H. Y. Tung, H. Strathmann, S. De, A. Ramdas, A. Smola, and A. Gretton. Generative models and model criticism via optimized maximum mean discrepancy. In International Conference on Learning Representations, 2017.
[71] M. F. Tipping and C. M. Bishop. Probabilistic principal component analysis. Journal of the Royal Statistical Society: Series B, 61 (3)):6111-63 Aug. 1999.
[72] D. Tran, M. Hoffman, D. Moore, C. Suter, S. Vasudevan, A. Radul, M. Johnson, and R. A. Saurous. Simple, distributed, and accelerated probabilistic programming. Tn Neural Information Processing Systems, 2018.
14

---

[73] A. van den Oord, S. Dieleman, H. Zen, K. Simonyan, O. Vinyals, A. Graves, N. Kalch- brenner, A. Senior, and K. Kavukeuoglu. WaveNet: A generative model for raw audio. Sept. 2016.
[74] L. van der Maaten and G. Hinton. Visualizing data using t-SNE. Journal of Machine Learning Research, 9:2579-2600 2603, 2008.
[75] A. W. van der Vaart. Asymptotic Statistics. 1998.
[76] R. Vershynin. High-Dimensional Probability: An Introduction with Applications in Data Science. 2020.
[77] P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser, J. Bright, S. J. van der Walt, M. Brett,
J. Wilson, K. J. Millman, N. Mayorov, A. R. J. Nelson, F. Jones,
Kern, Larson, C. J. Carey, T. Polat, Y. Feng, E. W. Moore, J. VanderPlas, D. Laxalde, J. Perktold, Cimrman, I. Henriksen, E. A. Quintero, C. R. Harris, A. M. Archibald, A. H. Ríbeiro, F. Pedregosa, P. van Mulbregt, and SciPy 1.0 Contributors. SciPy 1.0: Fundamental Algorithms for Scientific Computing Python. Nature Methods, 17:261 272, 2020.
[78] Y. Voichek and D. Weigel. Identifying genetic variants underlying phenotypic variation in plants without complete genomes. Nature Geneties, 52(5):534-540, 2020.
[79] E. N. Weinstein and D. S. Marks. A structured observation distribution for generative biological sequence prediction and forecasting. In Internat ional Conference on Machine Learning, 139, pages 11068 11079, PMLR, 2021.
15


Reviews:
[{"review_id": "x5AwDpnPKe2", "reviewer": "Reviewer_MWYx", "summary": "The authors propose a generative sequence model, Bayesian embedded autoregressive (BEAR), which is non-parametric as opposed to parametric autoregressive models to address the misspecification problem in neural networks. The BEAR model is a Bayesian Markov model with an embedded autoregressive model into its prior. Using parametric model output as a prior allows BEAR to be used for non-parametric density estimation which makes asymptotic consistency possible. They also show that the learned concentration hyperparameter for a prior dirichlet distribution represents the amount of misspecification in the embedded autoregressive model (or other neural network that was embedded) that can be used as a warning when the model can not be trusted. Then, they prove that absolute performance of the distribution predicted by some model can be measured through testing whether it matches the optimal BEAR model. \n\nThey theoretically prove that the BEAR model is more expressive than the naive AR model and that it converges to the true distribution as the size of the dataset increases. BEAR models are scalable and uncertainty-aware, which make them a useful generative model for full biological sequences. ", "questions": "", "limitations": "", "rating": 7, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "EDIT 2021-08-19: Changed score to 7.\n\nThe paper addresses the interesting problem of sequence generation with a convergence guarantee, and is scalable to whole genome sequences which make it practically used as a generative model. It is well-written and all claims are followed by theoretical proofs and simulation results. My main concern is that it does not reference or compare to the related literature on sequence compression.\n\nMajor comments:\n\n- What is the utility of a generative model of biological sequences? The authors say that it would be useful to have a generative model, but not what it might be used for. The closest experiment to being practical is in section \"Testing hypotheses\", in which the authors design a statistical test for whether two sets of sequences are generated from the same distribution. Is this the main aim of BEAR?\n- Density estimation is closely related to compression. There is a large literature of nucleotide sequence compression methods, such as the following, which the authors should cite and compare and contrast. \n\t- https://academic.oup.com/gigascience/article/9/11/giaa119/5974977\n\t- https://academic.oup.com/bib/article/16/1/1/240815\n- The manuscript includes a huge 62-page supplement. Much of the paper is unclear without reference to the supplement, including many of the experiments. The manuscript would be improved with a narrower focus, particularly on those tasks that are practically relevant. \n- In particular, the formulation of fk(θ) is not in the main manuscript, although it is very needed to know the parameters that are going to be learned during training. For example, what is B when β* > 0 in AR model if it is not learned during fitting? What is its value during sample generation? \n- In the “measuring misspecification” subsection in the results section, the relation between the misspecification and flexibility is not clear. The misspecification value has been shown to be proportional to a divergence between AR model and true distribution, however it is not clear how the divergence is related to the flexibility of the model.\n- In the two-sample test experiment (evaluating a simulation model, ART), it has been shown that differences between real and simulation data can be observed when the test is done for individual lags. The behavior of bayes factor in terms of lag is not clear, for example why it increases and then decreases when comparing simulation and real data.", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "Yes", "review_text": "", "extra_scores": {"ethical_concerns": "None", "needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "6"}, "extra_text": {}}, {"review_id": "s8hGSY3XNRS", "reviewer": "Reviewer_KMFd", "summary": "In this study, authors have proposed a Bayesian embedded autoregressive model to tackle the problem of measuring and making DNA. They have implemented the model using a parametric auto-regressive model. The model is programmed in PYTHON language. They try to overcome challenges of  having enough flexibility to account for genomic complexity, interpretability to reach scientific conclusions, and scalability to train on billions of nucleotides. The authors have explored, theoretically and empirically, applications of BEAR models to a variety of statistical problems  including density estimation, robust parameter estimation, goodness-of-fit tests, and two-sample tests. ", "questions": "", "limitations": "", "rating": 9, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "The paper is well-written in English, and impressively organized. The major novelty of the model comes from nonparametrically perturb a parametric model, studying its theoretical properties, and developing algorithms and implementations for terabyte-scale inference. The model outperforms other standard models, however, it has some limitations such as not explicitly account for paired-\nend read information, or other sources of long-distance information; But this model can be a starting point to address all these challenges. Authors have been honest about the limitations of their work. \n\n", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "Not applicable", "review_text": "", "extra_scores": {"ethical_concerns": "No concerns", "needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "6 hours"}, "extra_text": {}}, {"review_id": "UVDxyGDK4ie", "reviewer": "Reviewer_KhK3", "summary": "The paper introduces a generative Bayesian Markovian framework and a set of underlying models for DNA sequences to estimate the probability density of the current nucleotide given $L$ previous nucleotides $p(X_i|X_{i-L:i-1})$.\nThe authors investigate theoretically and empirically some applications of the framework such as density estimation of sequence data, parameter estimation and goodness-of-fit tests.", "questions": "", "limitations": "", "rating": 7, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "The paper is clearly written and well organized.\n\nThe novelty of the work lies in the theoretical description of the BEAR model and its applications to certain statistical problems for DNA modelling. The BEAR model is based on the Bayesian Markovian framework. It uses an embedded autoregressive (AR) model to estimate the hyperparamters of the Dirichlet prior distribution over the probability density of the kmers (nucleotide subsequences of length $L$). \nThe claims are adequately, theoretically and experimentally, supported.\n\n#### Comments\n\n1) I think the paper is too long for a conference venue. The submitted manuscript respects the 9 pages limit but each main section constantly refers the reader to other sections of the supplemental material. The presented models and the applications cannot be well understood without the support of this material. The theoretical results presented in the main text are a summary of a set of theorems and propositions about models, density estimation, parameter estimation and hypothesis testing that are detailed in 45 pages in the supplemental material. It is the same case for empirical results, which are supported by 17 pages in the supplemental material. \n\n2) It is not clear whether the finite-lag Markov models described in Section D are new or does this name refer to an already known family of Markov models? From their definition, they are quite similar to variable-order Markov models. If not, what is the difference between both of these models? The variable-order Markov models have been used in nucleotide sequence analysis, e.g., in virus genome subtyping [1]. \n\n3) How does this work (embedding a parametric model in a Dirichlet prior) relate to the work of Tzu-Tsung Wong on Generalized Dirichlet distribution in a Bayesian analysis [2, 3]?\n\n4) In the empirical evaluation, the authors designed “new” AR baseline models to embed them within the BEAR model and also to compare with. The results show that most of these baseline models are too weak and that their performance is lower or close to chance values. The comparison should also be done with other published methods or other deep architectures (RNNs, VAE) regardless of the AR framework (even RNNs could be considered as AR) [4, 5, 6]. In case the published methods cannot scale with terabyte-scale datasets (which I don't think is the case, e.g. metagenomics classification with linear methods [7] and UniRep for protein deep representation with mLSTM [8]), it will always be interesting to have comparison results on medium-size datasets. \n\n#### Minor comments\n\n1) Lines 23-25: Genome modelling using genetic variant matrix is usually used in GWAS studies. However, there is large repertoire of studies that models the whole genome at nucleotide resolution in different genomics analysis applications [4, 9] \n\n2) Line 60: Rephrase: Crucial to our the theoretical and empirical analysis is the statistical setting.\n\n#### References\n* [1] Struck, D., et al. COMET: adaptive context-based modeling for ultrafast HIV-1 subtype identification. 2014.\n* [2] Wong, T.T. Generalized Dirichlet distribution in Bayesian analysis. 1998.\n* [3] Wong, T.T. and Liu, C.R. An efficient parameter estimation method for generalized Dirichlet priors in naïve Bayesian classifiers with multinomial models. 2016.\n* [4] Im, Jinho, Byungkyu Park, and Kyungsook Han. A generative model for constructing nucleic acid sequences binding to a protein. 2019\n* [5] Liang, Qiaoxing, et al. DeepMicrobes: taxonomic classification for metagenomics with deep learning. 2020.\n* [6] Riesselman, Adam J., John B. Ingraham, and Debora S. Marks. Deep generative models of genetic variation capture the effects of mutations. 2018.\n* [7] Vervier, K., et al. Large-scale machine learning for metagenomics sequence classification. 2016.\n* [8] Alley, Ethan C., et al. \"Unified rational protein engineering with sequence-based deep representation learning. 2019.\n* [9] Lim, D., and Blanchette M. EvoLSTM: context-dependent models of sequence evolution using a sequence-to-sequence LSTM. 2020.\n", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "The limitations of the BEAR model are briefly discussed in Section 7.\nThe social impact of genetic-related models is highlighted in Section A.", "review_text": "", "extra_scores": {"ethical_concerns": "Not applicable", "needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "12"}, "extra_text": {}}, {"review_id": "8zTfqKi4S6l", "reviewer": "Reviewer_mk3X", "summary": "In this work, authors propose a new generative sequence model, called the Bayesian embedded autoregressive (BEAR) model. BEAR is a nonparametric Bayesian Markov model which embeds an autoregressive (AR) model into its prior. Theoretical guarantees are provided for statistical problems including density estimation, misspecification detection, and hypothesis testing, and experimental results indicate BEAR's superior performance on a variety of datasets.", "questions": "", "limitations": "", "rating": 6, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "\\textbf{Strength}: \nBEAR builds generative probabilistic models at the genome scale, which is able to detect and overcome misspecification problem. \\\\\nIn the paper, authors derive the asymptotic consistency guarantees and convergence rates for nonparametric density estimation, and prove consistency results for goodness-of-fit tests and two-sample tests, which provide theoretical guarantees of BEAR.\\\\\nIn the experimental evaluation, BEAR outperforms other autoregressive models on large-scale genome datasets.\n\n\\textbf{Weakness}: \nI have several concerns regarding to the proposed method. \nFirst, I wonder if BEAR can handle the long-range sequential data?\\\\ \nNext, I wonder if the parameters of BEAR are also interpretable? It would be great if the authors could discuss further about how to quantify the degree of misspecification introduced in the Section Robust parameter estimation, and how to select the optimal lag $L$. \\\\\nThird, will the model be able to handle prior knowledge about some special elements on the genome sequence, such as transposable elements? \n", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "The authors have discussed the limitations in the discussion section, and the proposed method does not have any negative societal impacts.", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "4"}, "extra_text": {}}]

Author rebuttal:

Meta-review: All four reviewers recommend accepting the submission, but three reviewers are not highly confident in their assessments and I am not highly confident in the other reviewer’s assessment; there is a disconnect between what this reviewer wrote and what rating this reviewer assigned. Therefore, I reviewed this submission carefully myself. The submission is well written, novel, and makes an important contribution to Bayesian modeling for genomic sequence data; I recommend accepting it.
Decision: Accept (Poster)