A Generative Nonparametric Bayesian Model for Whole Genomes
Year: 2021

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 o er a flexible and scalable framework, with theoretical guarantees, for building and critiquing generative models at the whole genome scale.

heading: 1 Introduction
text: Measuring and making DNA is central to modern biology and biomedicine. Generative probabilistic modeling o ers a framework for learning from sequencing data and forming experimentally testable predictions of unobserved or future sequences that can be synthesized in the laboratory [19, 31, 63]. Existing approaches to genome modeling typically preprocess the data to build a matrix of genetic variants such as single nucleotide polymorphisms [25, 57]. However, most modes of 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 [12, 44, 50, 69, 78]. Ideally, we would be able to directly model genome sequencing data and/or assembled genome sequences. However, building generative models that work with raw nucleotides, not matrices of alleles, raises the
úThese authors contributed equally.
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
. . .GTATCC A T G C
A T G C
A
T G
A
A T
G C AT ATATA
AG GC
AC CA
CT
TAC ACA CAC ACT CTA
TAG AGC GCA
B

heading: 2 Bayesian embedded autoregressive models
text: We first briefly review autoregressive (AR) models as applied to sequences of discrete characters. Let f(◊) denote an autoregressive function with parameter ◊ and let L denote the lag of the autoregressive model; then the AR model generates data as
Xi|Xi≠L:i≠1 ≥ Categorical(fXi≠L:i≠1(◊)), (1)
where i indexes position in the sequence X and Xi≠L:i≠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 a start symbol ÿ at the beginning and a stop symbol $ at the end of each sequence; letters Xi 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,
L ≥ fi(l), vk ≥ Dirichlet ! 1
h fk(◊)
"
for all k,
Xi|Xi≠L:i≠1 ≥ Categorical(vXi≠L:i≠1), (2)
where fi(l) is a prior on the lag with support up to infinity, h > 0 is a concentration hyperparameter, and k is a length L kmer. The BEAR model has three key properties (Fig. 1). First, the unrestricted transition parameter v and lag L allow the model to capture exact conditional distributions of pú to arbitrarily high order: pú(Xi|Xi≠1) at L = 1, then p
ú(Xi|Xi≠2, 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. 1). The optimal h provides a measurement of the amount of misspecification in the AR model (Section 4). Third, the choice of the conjugate Dirichlet prior allows the conditional marginals p((Xn)Nn=1|L, h, ◊) to be computed analytically, and (since L is one-dimensional) the total marginal likelihood p((Xn)Nn=1|h, ◊) to be estimated tractably. This allows BEAR models to be used for hypothesis testing (Section 5). There are 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 ◊. Let #(k, b) denote the number of times the length L kmer k is seen followed by the letter or stop symbol b in the dataset (Xn)Nn=1. Using a high-performance kmer counter optimized for nucleotide data, KMC, we can compute the count matrix #(·, ·) for all observed kmers k in terabyte-scale datasets, even when the matrix does not fit in main memory (Section J.2) [39]. To optimize h and ◊, we take advantage of the fact that the log conditional marginal likelihood can be written as a sum over observed kmers,
log p((Xn)Nn=1|L, h, ◊) = ÿ
k:#k>0 log
5 ( q b 1 h fkb(◊))
r b ( 1h fkb(◊))
r b ( 1h fkb(◊) + #(k, b)) ( q
b 1 h fkb(◊) + #(k, b))
6
. (3)
This decomposition lets us construct unbiased stochastic estimates of the gradient with respect to h and ◊ 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 https://github. com/debbiemarkslab/BEAR. 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 fk(◊) depends on k linearly as a function of both individual positions and pairwise interactions between positions, with the strength of the pairwise interaction weighted by a parameter —ú (Section I.1). We first fit (using maximum likelihood) a linear AR model that lacks pairwise terms and is thus misspecified when —ú > 0. Since the AR model is misspecified, it does not asymptotically approach the true data-generating distribution pú (Fig. 2A, gray). We next computed the posterior of a vanilla BEAR model without the embedded AR in its prior, instead using the Je reys prior vk ≥iid Dirichlet(1/2, . . . , 1/2). The vanilla BEAR model asymptotically 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 a BEAR model with the misspecified linear AR model embedded, using our empirical Bayes procedure. The BEAR model performs just as 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. 2A, blue and yellow). When the AR model is well-specified, the empirical Bayes estimates of the parameters ◊ under the BEAR model match the maximum likelihood estimates of ◊ under the AR model
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. 2B). This warning emerges early: h 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 parametric model [48], following a similar strategy to the Polya tree method proposed by Berger and Guglielmi [6]. As in Berger and Guglielmi [6], we use Dirichlet priors centered at the parametric model’s predictions, and construct tractable goodness-of-fit tests by exploiting Dirichlet-categorical conjugacy. BEAR models extend these ideas from one-dimensional continuous data to finite-length sequences of discrete characters. Markov and AR models have a long history and wide range of applications in biological sequence analysis [21, 51, 58]. Compression methods, in particular, often rely on accurate density estimation and use Markov or AR models to achieve it [17, 55, 56, 67]. We establish theoretical guarantees for density estimation with fully Bayesian Markov 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 a sequence dataset; in the limit h æ Œ, 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 [11, 33]. Comparisons between genomes and other sequences are often made on the basis of kmer counts; our two-sample test provides a generative perspective on this idea [3, 16, 78]. BEAR models are also connected to ideas in natural language processing, where kmers are referred to as ngrams. Under the vanilla BEAR model, the mean of the posterior predictive distribution conditional on L corresponds to an ngram additive smoothing model [9]. Comparisons between datasets using their ngram counts are also common in model evaluation metrics such as the BLEU score [53].

heading: 3 Density estimation
text: The density estimation problem is that of estimating pú given data (Xn)Nn=1 drawn i.i.d. from pú. Density estimation is particularly crucial for biological sequence analysis due to its connections to fitness estimation [31, 65]. State-of-the-art mutation e ect prediction methods and clinical variant interpretation methods rely on density estimates of evolutionary sequence data [20, 60]. Density estimation with generative models is particularly useful for protein design, as samples from accurate density estimates are likely to be functional and
can be synthesized in the laboratory [63, 66]. 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 [79]. Here, we show that the posterior distribution of the BEAR model is consistent and will concentrate on pú as N æ Œ, regardless of what pú actually is, so long as pú generates finite length sequences almost surely (a.s.). We first study the expressiveness of BEAR models. Let ML be the set of Markov models pv with transition probabilities v and lag L that generate finite length strings a.s.. Note that M1 µ M2 µ . . .. Define the union M = fiŒL=1ML. 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 1-4 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 data-generating distributions, they can approximate pú arbitrarily well as L 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. [23]. The most important assumption is that p
ú is subexponential, meaning that for some t > 0, Epú exp(t|X|) < Œ where |X| is the sequence length. Let (·|(Xn)Nn=1) denote the posterior over sequence distributions. Let B(pú, ”) denote a ball of radius ” centered at pú, using the Hellinger distance. Summary of Theorem 35 Given M > 0 large enough and ‘ œ (0, 1) small enough, we have (B(pú, MN≠ 12 ‘)|(Xn)Nn=1) æ 1 in probability. A proof is in Section H and simulations in Section I.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 a rate of convergence: in a parametric model, the uncertainty would shrink as N≠ 12 , but here the rate is slower, N≠ 12 ‘, a price paid for the nonparametric model’s expressivity [23, 29, 35]. The proof includes a variety of new theoretical constructions and algorithms that are used to approximate subexponential sequence distributions.

heading: 4 Robust parameter estimation
text: To derive a biological understanding of mutational processes, evolutionary history, functional 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 [34]. To reach robust scientific conclusions, therefore, parameter estimates should ideally come with a 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 ◊, as well as the hyperparameter h, showing that h diagnoses misspecification in the embedded AR model. Our analysis builds o the study of empirical Bayes consistency in Petrone et al. [54], 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 arbitrary and large value [30]. Define pú(L) = argminpvœMLkl(púÎpv) as the closest model in ML to p
ú, and define vú such that pvú = pú(L) (note pú(L) æ pú as L æ Œ). We say that the AR model is misspecified “at resolution L” if f cannot approximate pú(L), i.e. if there does not exist some sequence of parameter values ◊̃N such that pf(◊̃N ) æ p
ú(L) as N æ Œ; otherwise, the AR model is well-specified at resolution L. Now we can study empirical Bayes estimates of h and ◊, denoted hN and ◊N . Summary of Propositions 15-20 Let (hN )ŒN=1 and (◊N )ŒN=1 be sequences maximizing the BEAR marginal likelihood p((Xn)Nn=1|L, h, ◊) for each N . If the model is well-specified at resolution L, then hN N1/4≠‘ æ 0 for every ‘ > 0 and pf(◊N ) æ pú(L) in distribution, with both sequences converging in probability. On the other hand, if the model is misspecified at
resolution L, then hN is eventually bounded below by some positive (non-zero) number a.s.. Proofs are in Section F and simulations in Section I.1. The implication of this result is that when the AR model is well-specified, hN converges to zero (at a rate that is a power of the dataset size) and ◊N converges to the parameter value ◊ú at which the AR model matches the data (Corollary 16). On the other hand, when the AR model is misspecified, hN does not converge to zero; heuristically, we find instead that hN is approximately proportional to a divergence between pú(L) and the AR model,
hN Ã ÿ
kœaccL(pú)
Q
akl(fk(◊N )Îvúk) + log(N) ÿ
b/œsuppL(pú)|k
fk,b(◊N )
R
b , (4)
where accL(pú) = {k | pú(#k > 0) > 0} is the set of kmers with non-zero probability and suppL(pú)|k = {b | pú(#(k, b) > 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 a parameter estimate ◊N , a value hN which tells you the amount (from zero to infinity) of misspecification in the AR model. If hN is close to zero, you can trust the estimate ◊N .

heading: 5 Hypothesis testing
text: 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 [45]. 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 statistics, such as average amino acid hydrophobicity or sequence length [63, 66]. Given a dataset (Xn)Nn=1 ≥ pú i.i.d., a goodness-of-fit test asks whether or not the data distribution pú matches a model distribution p̃. It takes into account all possible distributions pú including those that di er from p̃ in a manner that cannot be captured by finitely many summary statistics. We propose a goodness-of-fit test that compares the null hypothesis H0 : pú = p̃ to the alternative H1 : pú ”= p̃ using the Bayes factor BF = p((Xn)Nn=1|h, ◊)/p̃(X1:n), where p((Xn)Nn=1|h, ◊) = q
L p((Xn)Nn=1|L, h, ◊)fi(L) is the marginal likelihood under the BEAR model. Note that practically, the sum over L is straightforward to approximate by truncation, and that the test can be computed in time linear in the amount of data. 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 fi(L) = 0 for L larger than some arbitrary L̃ but fi(L) > 0 for L Æ L̃ [30]. We treat ◊ and h > 0 as fixed. Summary of Proposition 21 If p̃ is at least as close to pú as pú(L) is, as measured by kl(púÎ·), then BF æ 0 in probability as N æ Œ. On the other hand, if pú(L) is closer than p̃, then BF æ Œ in probability. A proof is in Section G.1 and simulations in Section I.3. An important practical limitation on nonparametric 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 BF(L, h, ◊) = p((Xn)Nn=1|L, h, ◊)/p̃((Xn)Nn=1) with any choice of L, h > 0, and ◊. Thus in practice to increase power we can maximize the value of BF(L, h, ◊) as a function of L, h, and/or ◊ (note that this approach is heuristic, since we have not proven the consistency of the maximized Bayes factor). Berger and Guglielmi [6] provide extensive methodological guidance on using analogous tests constructed with Polya trees. Based on their recommendations, we suggest first choosing ◊ such that pf(◊) is as close as possible to p̃, then plotting the Bayes factor as a function of h and/or L to identify the maximum value and confirm that any conclusion is robust to changes in h and/or L. Another challenge in nonparametric hypothesis testing is that it can be di cult 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 [43, 70].
Two-sample test A two-sample test asks whether or not two datasets (Xn)Nn=1 and (X Õn)N Õ
n=1 are drawn from the same distribution. E orts to compare di erent sequence datasets are widespread in biology: for instance, researchers often wish to determine whether two microbiome samples, taken under di erent conditions or at di erent timepoints, are the same up to sampling noise [44]. Two-sample tests can also be used to evaluate generative sequence models that lack tractable likelihoods (for which the goodness-of-fit test proposed above does not apply) such as energy-based models or implicit models like GANs and biophysical simulators [27, 42, 49]. Assume (Xn)Nn=1 ≥ p1 and (X Õn)N Õ
n=1 ≥ p2 i.i.d.. Our BEAR test compares the null hypothesis H0 : p1 = p2 to the alternative H1 : p1 ”= p2 using the Bayes factor BF = p((Xn)Nn=1|h, ◊)p((X Õn)N Õ n=1|h, ◊)/p((Xn)Nn=1, (X Õn)N Õ
n=1|h, ◊). As in the goodness-of-fit 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 ◊. Summary of Proposition 22 If p(L̃)1 = p (L̃) 2 , then BF æ 0 as N æ Œ in probability. Otherwise, if p(L̃)1 ”= p (L̃) 2 , then BF æ Œ in probability. A proof is in Section G.2 and simulations in Section I.3.

heading: 6 Results
text: Predicting sequences We sought to evaluate BEAR models as compared to AR models on the task of predicting real nucleotide (nt) sequences. We considered eleven datasets of four di erent 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 ≠ 1010 nt and in individual sequence length from ≥ 102 ≠ 106 nt (Table S1). 25% of data was randomly held out for testing, in the form of entire sequences (reads, genomes, etc., see Table S2); 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 a linear AR model and a deep convolutional neural network (CNN) AR model with > 10◊ more parameters, both of which are common models used across a range of applications; we also designed a biologically-structured AR model which makes predictions based on a reference genome and a Jukes-Cantor mutation model
(Section L.1) [56, 67]. 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 1) and accuracy (Table S3) in all datasets, even when the model lag L is held fixed for comparison (Section L.3). 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. S13). Remarkably, the optimal lag L chosen by empirical Bayes is often quite short, less than 20 nt (Table S4). The improvements o ered 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 3◊ instead of 100◊, the improvement in prediction accuracy would have been greater than 10 percentage points instead of 0.1 (Section L.4; Fig. S14).
Table 2: Diagnostic h. Abbreviations as in Table 1.
Dataset Lin. CNN Ref. YSD1 5.528 5.461 4.183 A. th. 1 2.765 2.756 2.990 A. th. 2 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 -
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. Examining 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 h 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 su ers 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 scientific models outperform deep neural networks it is thanks to their low variance in the small data regime.
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 an alternative non-probabilistic extrapolation method that is widely used in biology, local assembly (Fig. 3A; Section M). In this example the assembly algorithm SPAdes returns four possible assemblies, a relatively large number compared to other reads in the dataset (Fig. 3A stars) [5]. Samples from the BEAR model include these four possibilities, but also many more, some with higher probability. The distribution over possible nucleotide choices under the BEAR model is much wider than the number of assemblies would suggest: it has a 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. 3B). These observations suggest that SPAdes, which does not provide a measurement of uncertainty, may not be capturing the full range of possible sequences. Visualizing data Methods for learning local representations or features of biological sequences can be powerful tools for visualization and semisupervised learning [7]. One approach to extracting such representations is to learn a generative model q(X1, . . . , XL+1) of kmers, for instance using a variational autoencoder. While such models are not autoregressive, the small size of the DNA alphabet makes it tractable to estimate the conditional q(XL+1|X1:L) by Bayes’ rule, and this conditional can then be embedded into a BEAR model. We applied this strategy to probabilistic 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. 3CD; 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.
Testing hypotheses The question of when and how microbiomes change is widespread, but has in the past relied on summary statistics of sequencing datasets [44]. Schreiber et al. [64] studied changes in patient urine microbiomes before and after kidney transplant, and performed both unbiased metagenomic sequencing and diagnostic quantitative polymerase chain reaction (qPCR) for a specific virus associated with complications (JC polyomavirus). They found evidence of donor-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 Je reys 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 S6; Section O.1). These results show, in a small example, that the two-sample test has su cient power to detect microbiome changes in real data, and can be consistent with more specific tests. We next applied BEAR hypothesis tests to 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 > 108) for misspecification in each A. thaliana whole genome sequencing dataset, and this conclusion is robust to a wide range of h values (Fig. 3E; Section O.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 di erences (Fig. 3F; Section O.2). Note that as L increases, there is a tradeo : tests with larger lag can detect more subtle di erences between the two distributions, but have less statistical power since they must consider a larger set of possible distributions. Thus the
Bayes factor first increases and then decreases with lag, reaching a peak at intermediate values where there is the most evidence of di erence. To understand in detail the source of the detected di erences 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. 3G). One group mapped to chloroplasts, an organelle 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 a 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.

heading: 7 Discussion
text: In this article we proposed the nonparametric 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-probabilistic 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 pairedend 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.

Review summaries:
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. 

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.   They 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. 

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.