# Nonparametric Bayes Differential Analysis of Multigroup DNA Methylation Data

Chiyu Gu
Bayer Crop Science, Chesterfield, Missouri
and
Veerabhadran Baladandayuthapani
Department of Biostatistics, University of Michigan
and
Subharup Guha
Department of Biostatistics, University of Florida

###### Abstract

DNA methylation datasets in cancer studies are comprised of measurements on a
large number of genomic locations called cytosine-phosphate-guanine (CpG)
sites with complex correlation structures. A fundamental goal of these studies
is the development of statistical techniques that can identify disease genomic
signatures across multiple patient groups defined by different experimental or
biological conditions. We propose BayesDiff, a nonparametric Bayesian approach
for differential analysis relying on a novel class of first order mixture
models called the Sticky Pitman-Yor process or two-restaurant two-cuisine
franchise (2R2CF). The BayesDiff methodology flexibly utilizes information
from all CpG sites or probes, adaptively accommodates any serial dependence
due to the widely varying inter-probe distances, and performs simultaneous
inferences about the differential genomic signature of the patient groups.
Using simulation studies, we demonstrate the effectiveness of the BayesDiff
procedure relative to existing statistical techniques for differential DNA
methylation. The methodology is applied to analyze a gastrointestinal (GI)
cancer dataset that displays both serial correlations and interaction
patterns. The results support and complement known aspects of DNA methylation
and gene association in upper GI cancers.

Keywords: 2R2CF; First order models; Mixture models; Sticky Pitman-Yor
process; Two-restaurant two-cuisine franchise

##  1 Introduction

Recent advances in array-based and next-generation sequencing (NGS)
technologies have revolutionized biomedical research, especially in cancer.
The rapid decline in the cost of genome technologies has facilitated the
availability of datasets involving intrinsically different sizes and scales of
high-throughput data and provided genome-wide, high resolution information
about the biology of cancer. A common analytical goal is the identification of
differential genomic signatures between groups of samples corresponding to
different treatments or biological conditions, e.g., treatment arms, response
to adjuvant chemotherapy, tumor subtypes, or cancer stages. The challenges
include the high dimensionality of genomic biomarkers or probes, usually in
the hundreds of thousands, and the relatively small number of patient samples,
usually no more than a few hundred. This “small $n$, large $p$” setting
results in unstable inferences due to collinearity. Further, there exist
complex interaction patterns, such as signaling or functional pathway-based
interactions, and location-based serial correlation for high-throughput
sequencing data. These data attributes significantly affect the reliability of
statistical techniques for detecting differential genomic signatures.

##### Differential DNA methylation in cancer studies

DNA methylation is an important epigenetic mechanism that involves the
addition of a methyl (CH3) group to DNA, resulting in the modification of gene
functions. It typically occurs at specific genomic locations called cytosine-
phosphate-guanine (CpG) sites. Alterations in DNA methylation, e.g.,
hypomethylation of oncogenes and hypermethylation of tumor suppressor genes,
are often associated with the development and progression of cancer (Feinberg
& Tycko, 2004). It was previously believed that these alterations occur almost
exclusively at promoter regions known as CpG islands, i.e., chromosomal
regions with high concentrations of CpG sites. However, with the advent of
high-throughput technologies, it has been shown that a significant proportion
of cancer-related alterations do not occur in promoters or CpG islands
(Irizarry et al., 2009), prompting higher resolution, epigenome-wide
investigations.

Gastrointestinal (GI) cancer, the most common form of cancer in the U.S.
(Siegel et al., 2017), refers to malignant conditions affecting the digestive
system associated with epigenetic alterations (Vedeld et al., 2017). Molecular
characterization of different cancer types, facilitated by the identification
of differentially methylated CpG sites, is therefore key to better
understanding GI cancer. In the motivating application, we analyze methylation
profiles publicly available from The Cancer Genome Atlas (TCGA) project,
consisting of 1,224 tumor samples belonging to four GI cancers of the upper
digestive tract: stomach adenocarcinoma (STAD), liver hepatocellular carcinoma
(LIHC), esophageal carcinoma (ESCA) and pancreatic adenocarcinoma (PAAD). For
485,577 probes, where each probe is mapped to a CpG site, DNA methylation
levels or Beta-values ranging from 0 (no methylation) to 1 (full methylation)
were measured using the Illumina Human Methylation 450 platform.

Figure 1: Methylation levels of CpG sites near gene TP53 for randomly chosen
tumor samples of the TCGA upper GI dataset. Each plotted point represents the
methylation level at a CpG site, with the shapes and colors corresponding to
different GI cancers indicated in the legend. The vertical dashed lines
demarcate the TP53 gene boundaries.

Figure 1 displays the methylation levels for CpG sites near TP53, a tumor
suppressor gene located on chromosome $17$. A random subset of the tumor
samples was chosen to facilitate an informal visual evaluation. Each plotted
point represents the methylation level of a tumor sample at a CpG site. As
indicated in the figure legend, the four sets of colors and shapes of the
points represent the four upper GI cancers. The vertical dashed lines indicate
the boundaries of the TP53 gene. Although differential methylation is clearly
visible at some CpG sites, the differences are generally subtle, demonstrating
the need for sophisticated statistical analyses. An obvious feature is the
correlation between the apparent methylation statuses of nearby CpG sites
(Eckhardt et al., 2006; Irizarry et al., 2008; Leek et al., 2010). The
dependence of proximal CpG sites is also seen in Figure 6 of Supplementary
Material, where we find high first order autocorrelations and highly
significant tests for serial correlations. Furthermore, the variability of the
inter-probe spacings in Figure 1 suggests that the need to model distance-
based dependencies.

##### Existing statistical approaches for differential DNA methylation and
limitations

Numerous frequentist and Bayesian methods have been developed for differential
DNA methylation, and can be broadly classified into four categories: (i)
Testing-based methods, such as Illumina Methylation Analyzer (IMA) (Wang et
al., 2012), City of Hope CpG Island Analysis Pipeline (COHCAP) (Warden et al.,
2013), and BSmooth (Hansen et al., 2012). These methods rely on two-sample or
multiple-sample tests for the mean group differences at each CpG site. (ii)
Regression based models, such as Methylkit (Akalin et al., 2012), bump hunting
(Jaffe et al., 2012), Biseq (Hebestreit et al., 2013), and RADMeth (Dolzhenko
& Smith, 2014). After applying smoothing or other adjustments, these methods
fit individual regression models for each CpG site and test for significance.
(iii) Beta-binomial model-based methods, such as MOABS (Sun et al., 2014), DSS
(Feng et al., 2014), and methylSig (Park et al., 2014). These methods fit
separate models to each CpG site. (iv) Hidden Markov models (HMMs), such as
MethPipe (Song et al., 2013), Bisulfighter (Saito et al., 2014), and HMM-DM
(Yu & Sun, 2016). These methods detect differentially methylated sites based
on inferred hidden states.

The aforementioned methods have several deficiencies. By fitting separate
models to each probe, most methods ignore the strong correlations between
neighboring probes. This reduces the detection power for the relatively small
sample sizes. Additionally, beta-binomial, HMM, and most testing-based methods
are able to accommodate only two treatments and rely on inefficient
adjustments to compare multiple treatments. The methods that account for
serial dependence (e.g., HMMs) do not adjust for the widely varying distances
between the probes, and instead, assume uniform inter-probe dependencies. The
few methods that account for inter-probe distances (e.g., Hansen et al., 2012;
Jaffe et al., 2012; Hebestreit et al., 2013) rely on ad hoc parameter-tuning
procedures that do not adjust for distinctive data characteristics.

Motivated by these challenges, we propose general and flexible methodology for
differential analysis in DNA methylation data, referred to as BayesDiff.
Rather than fitting a separate model to each CpG site or probe, BayesDiff
relies on a global analytical framework for simultaneous inferences on the
probes that adapts to the unique data attributes. To diminish collinearity
effects and achieve dimension reduction, the probes are allocated to a
smaller, unknown number of latent clusters based on the similarities of
probes-specific multivariate parameters. Finally, differential state variables
of the probes delineate the genomic signature of the disease to fulfil the
main inferential goal.

For realistically modeling the probe-cluster allocation mechanism of DNA
methylation profiles, we devise an extension of Pitman-Yor processes (PYPs)
(Perman et al., 1992) called the Sticky PYP (equivalently, the two-restaurant
two-cuisine franchise). In addition to accounting for long-range biological
interactions, this nonparametric process accommodates distance-based serial
dependencies of the probes. Separately for the differential and non-
differential probes, it flexibly permits the data to direct the choice between
PYPs, and their special case, Dirichlet processes, in finding the best-fitting
allocation schemes.

We implement an inferential procedure for Sticky PYPs using a Markov chain
Monte Carlo (MCMC) algorithm specifically designed for posterior inferences in
the typically large methylation datasets. Simulation results show that our
approach significantly outperforms existing methods for multigroup comparisons
in data with or without serial correlation. For the motivating TCGA dataset,
in addition to confirming known features of DNA methylation and disease-gene
associations, the analysis reveals interesting aspects of the biological
mechanisms of upper GI cancers.

The rest of the paper is organized as follows. Section 2 describes the
BayesDiff approach, with Section 2.1 introducing the Sticky PYP or two-
restaurant two-cuisine franchise (2R2CF) for differential DNA methylation.
Section 3 outlines an effective inference procedure for detecting differential
probes. Section 4 uses artificial datasets with varying noise and correlation
levels to assess the accuracy of BayesDiff in detecting disease genomic
signatures and makes comparisons with established techniques for DNA
methylation data. The motivating upper GI dataset is analyzed using the
BayesDiff procedure in Section 5. Finally, conclusions and future related work
are discussed in Section 6.

##  2 The BayesDiff Model

Sequencing technologies measure DNA methylation levels of $p$ biomarkers
represented by CpG sites (“probes”) and $n$ matched patient or tissue samples
(“individuals”). Usually, $p$ is much larger than $n$. The methylation levels,
which belong to the interval $[0,1]$, are arranged in an $n\times p$ matrix of
proportions, $\mathbf{X}=((x_{ij}))$, for individuals $i$ and probes $j$, with
the probes sequentially indexed by their genomic locations. The distances
between adjacent probes are denoted by $e_{1},\ldots,e_{p-1}$, and typically
exhibit high variability. For instance, in the upper GI TCGA dataset, the
inter-probe distances range from 2 base pairs to a million base pairs; a base
pair is a unit of DNA length consisting of two nucleobases bound to each other
by hydrogen bonds (e.g., Baker et al., 2008).

Each individual $i$ is associated with a known experimental or biological
condition (“treatment”) denoted by $t_{i}$ and taking values in
$\\{1,...,T\\}$ with $T\geq 2$. In the motivating TCGA data, there are $T=4$
upper GI cancer types. We model the logit transformation of the methylation
levels, $z_{ij}=\log\left(x_{ij}/\left(1-x_{ij}\right)\right)$, as follows:

$z_{ij}\sim N\left(\xi_{i}+\chi_{j}+\theta_{t_{i}j},\sigma^{2}\right)$ (1)

where $\xi_{i}$ represents the $i$th subject’s random effect, $\chi_{j}$
represents the $j$th probe’s random effect, and $\theta_{tj}$ is the random
treatment $t$–probe $j$ interaction effect. The logit-transformed methylation
levels differ from M-values (Du et al., 2010), commonly used in differential
analyses, by $1-\log(2)=0.306$; however, the key results are identical.

The main inferential goal is the detection of differential probes, i.e.,
probes $j$ for which the elements of column vector $\boldsymbol{\theta}_{j}=$
$\left(\theta_{1j},\dots,\theta_{Tj}\right)^{\prime}$ are not all identical.
Consequently, we define a binary differential state variable, $s_{j}$, with
$s_{j}=1$ indicating that probe $j$ is not differential and $s_{j}=2$
indicating that it is differential:

$\displaystyle s_{j}=\begin{cases}1\qquad\text{if
}\theta_{1j}=\theta_{2j}=\dots=\theta_{Tj},\\\
2\qquad\text{otherwise},\end{cases}$ (2)

for $j=1,\dots,p$. Thus, the key parameters are $s_{1},\ldots,s_{p}$, with the
differential genomic signature consisting of the probes for which state
$s_{j}=2$. Motivated by the distance-dependent correlations of DNA methylation
data and deficiencies of existing statistical approaches, this paper fosters a
Bayesian nonparametric framework for random effects
$\boldsymbol{\theta}_{1},\ldots,\boldsymbol{\theta}_{p}$ determining the
differential state variables.

### Modeling probe clusters

As mentioned, in addition to high-dimensionality, the analytical challenges
include pervasive collinearity due to dependencies between physically proximal
probes. Additionally, there are long-range dependencies between non-adjacent
probes because of biological interactions, e.g., signaling or functional
pathways. To accommodate these dependence structures and extract information
from the large number of probes, we allocate the $p$ probes to a much smaller
number, $q$, of latent clusters based on the similarities of the random
effects $\boldsymbol{\theta}_{j}$. We favor clustering to dimension reduction
methods such as principal components analysis (PCA). Because each principal
component is a linear combination of all $p$ biomarkers, PCA is less useful in
cancer research because of its inability to select features, i.e., probes. By
contrast, our approach facilities biological interpretations by identifying
CpG sites relevant to the differential genomic signatures between the
treatments.

Suppose that an allocation variable, $c_{j}$, assigns probe $j$ to one of $q$
latent clusters, where $q$ is unknown. The event $[c_{j}=k]$ signifies that
the $j^{th}$ probe belongs to the $k^{th}$ latent cluster, $k=1,\dots,q$. We
assume that the $q$ clusters are associated with latent vectors,
$\boldsymbol{\lambda}_{1},\dots,\boldsymbol{\lambda}_{q}$, where the probe-
specific random effects and cluster-specific latent vectors have the relation:

$\boldsymbol{\theta}_{j}=\boldsymbol{\lambda}_{k}\quad\text{if $c_{j}=k$.}$
(3)

That is, all probes within a cluster are assumed to have identical random
effects equal to that cluster’s latent vector. The differential state
variables of the probes, defined in equation (2), then become a shared
attribute of their parent cluster, and clusters as a whole are either
differentially or non-differentially methylated. Further, if probe $j$ belongs
to cluster $k$ (i.e., $c_{j}=k$), then the condition
$\theta_{1j}=\theta_{2j}=\dots=\theta_{Tj}$ in equation (2) is equivalent to
$\lambda_{1k}=\lambda_{2k}=\dots=\lambda_{Tk}$, and differential clusters
constitute the set

$\mathcal{D}=\left\\{k:\lambda_{tk}\neq\lambda_{t^{{}^{\prime}}k},\text{ for
some }t\neq t^{{}^{\prime}},k=1,\ldots,q\right\\}.$ (4)

##### Mixture models for allocation

Bayesian infinite mixture models are a natural choice for allocating $p$
probes to a smaller, unknown number of latent clusters based on their random
effects similarities. Dirichlet processes (Ferguson, 1973a) are arguably the
most frequently used infinite mixture models; see Müller & Mitra (2013, chap.
4) for a comprehensive review. The use of Dirichlet processes to achieve
dimension reduction has precedence in the literature, albeit in unrelated
applications (see Medvedovic et al., 2004; Kim et al., 2006; Dunson et al.,
2008; Dunson & Park, 2008; Guha & Baladandayuthapani, 2016). Lijoi, Mena &
Prünster (2007a) advocated the use of Gibbs-type priors (Gnedin & Pitman,
2006) for accommodating more flexible clustering mechanisms and demonstrated
the utility of Pitman-Yor processes (PYPs) in genomic applications. An
overview of Gibbs-type priors and characterization of the learning mechanism
is provided by De Blasi et al. (2013). Formally, the PYP (Perman et al., 1992)
relies on a discount parameter $d\in[0,1)$, positive mass parameter $\alpha$,
and $T$-variate base distribution $W$, and is denoted by
$\mathcal{W}(d,\alpha,W)$. The value $d=0$ yields a Dirichlet process with
mass parameter $\alpha$ and base distribution $W$. Suppose
$\boldsymbol{\theta}_{1},\ldots,\boldsymbol{\theta}_{p}$ are distributed as
$\mathcal{W}(d,\alpha,W)$. The stick-breaking representation of
$\mathcal{W}(d,\alpha,W)$ is $\boldsymbol{\theta}_{j}\stackrel{{\scriptstyle
iid}}{{\sim}}\mathcal{P}$, where random distribution $\mathcal{P}$ is the
discrete mixture $\sum_{v=1}^{\infty}\pi_{v}\delta_{\boldsymbol{\phi}_{v}}$,
with $\delta_{\boldsymbol{\phi}_{v}}$ denoting a point mass located at the
atom $\boldsymbol{\phi}_{v}\stackrel{{\scriptstyle iid}}{{\sim}}W$. The random
stick-breaking probabilities have the form $\pi_{1}=V_{1}$, and
$\pi_{h}=V_{h}\prod_{v=1}^{h-1}(1-V_{v})$ for $h>1$, where
$V_{h}\stackrel{{\scriptstyle indep}}{{\sim}}\text{beta}(1-d,\alpha+hd)$. Guha
& Baladandayuthapani (2016) introduced VariScan, a technique that utilizes
PYPs and Dirichlet processes for clustering, variable selection, and
prediction in high-dimensional regression problems in general, and in gene
expression datasets in particular. They also demonstrated that PYPs are
overwhelmingly favored to Dirichlet processes in gene expression datasets,
which typically exhibit no serial correlation.

##### Limitations of existing mixture models

Although the aforementioned mixture models achieve dimension reduction and
account for the long-range biological interactions between non-adjacent
probes, a potential drawback is their implicit assumption of a priori
exchangeability of the probes. Consequently, these techniques cannot account
for the serial correlation in methylation data. Infinite HMMs, such as the
hierarchical Dirichlet process hidden Markov model (HDP-HMM) (Teh et al.,
2006) and Sticky HDP-HMM (Fox et al., 2011), may be utilized to fill this gap.
Although these models are a step in the right direction, they have several
undesirable features for differential analysis. First, the degree of first
order dependence is uniform irrespective of the inter-probe distances. This is
unrealistic in methylation datasets where the correlation typically decreases
with inter-probe distance (Hansen et al., 2012; Jaffe et al., 2012; Hebestreit
et al., 2013). Second, an ad hoc exploratory analysis of the GI cancer dataset
reveals that the serial correlation in the treatment-probe effects is weaker
than the serial dependence between the differential state variables defined in
equation (2). Although there may not be a biological explanation for this
phenomenon, this makes sense from a statistical perspective because the
differential states are binary functions of the treatment-probe interactions;
the differential states are more sensitive in detecting first order dependence
even when the higher-dimensional (and noisier) treatment-probe interactions
show little or no correlation. This suggests that a hypothetical two-group
Markov model, rather than an infinite-group Markov model such as HDP-HMM or
Sticky HDP-HMM, would provide a better fit for the data. Third, the range of
allocation patterns supported by infinite HMMs is relatively limited. In
particular, realistic allocation patterns, such as power law decays in the
cluster sizes and large numbers of small-sized clusters, a common feature of
cancer datasets (Lijoi et al., 2007b), are assigned relatively small prior
probabilities by infinite HMMs.

###  2.1 Sticky PYP: A Two-restaurant, Two-cuisine Franchise (2R2CF) for
Differential Analysis

We invent a mixture model called the Sticky PYP. Essentially, this is a cohort
of regular PYPs that generates the probe-specific random effects by switching
generative PYPs at random locations along the probe sequence. Alternatively,
the well-known Chinese restaurant franchise (CRF) metaphor for HDP-HMMs and
Sticky HDP-HMMs can be generalized to the two-restaurant two-cuisine franchise
(2R2CF) to provide an equivalent representation of Sticky PYPs appropriate for
differential analysis. We first present a descriptive overview of 2R2CF.

Imagine that a franchise has two restaurants, labeled 1 and 2. Each restaurant
consists of two sections, also labeled 1 and 2. Each section serves a single
cuisine and its section-cuisine menu consists of infinite dishes. Section 1 at
both restaurants exclusively serves cuisine 1. The cuisine 1 menus, along with
selection probabilities associated with the dishes, are identical at the two
restaurants. Similarly, section 2 at both restaurants exclusively serves
cuisine 2, and the cuisine 2 menus are also identical at the restaurants.

A succession of $p$ customers, representing the CpG sites or probes, arrive at
the franchise. The waiting times between successive customers represent the
inter-probe distances, $e_{1},\dots,e_{p-1}$. Each customer selects a
restaurant and then selects a section (equivalently, cuisine) in the
restaurant. Each restaurant section has an infinite number of tables, and a
customer may either sit at a table already occupied by previous customers, or
sit at a new table. All customers at a table are served the same dish,
previously chosen from the section-cuisine menu by the first customer who sat
at that table. A customer who sits at a new table is allowed to independently
pick a dish from the cuisine menu with a cuisine-specific probability
associated with each dish. As a result, multiple tables in a restaurant
section may serve the same dish.

Restaurant 1 specializes in cuisine 1. Consequently, section 1 is more popular
with restaurant 1 patrons. Similarly, restaurant 2 specializes in cuisine 2,
and so, restaurant 2 customers tend to favor section 2 over section 1. By
design, if a customer has eaten a cuisine 1 (2) dish, then the next customer
is more likely to visit restaurant 1 (2), where cuisine 1 (2) is more popular.
In this manner, each customer tends to select the same cuisine as the previous
customer.

In the metaphor, cuisine 1 symbolizes the non-differential state and cuisine 2
symbolizes the differential state. The dish that franchise customer $j$ eats
represents the probe-specific random effect, $\boldsymbol{\theta}_{j}$. Since
cuisine 1 represents the non-differential state, its dishes are characterized
by $T$-variate random vectors with all equal elements; see equation (2). In
contrast, cuisine 2 (differential state) dishes are characterized by
$T$-variate random vectors with at least two unequal elements.

The dependence in the restaurant and cuisine choices of consecutive customers
account for long runs of differential or non-differential states in DNA
methylation data. However, a customer’s influence on the next customer
diminishes as the time interval between the two customers increases. That is,
the differential statuses of two adjacent probes become statistically
independent in the limit as the inter-probe distance grows.

The 2R2CF process is illustrated in Figure 2 and discussed below in greater
detail. The following specification conditions on $G$, an unknown distribution
in $\mathcal{R}$ that is assigned a Dirichlet process prior with mass
parameter $\beta>0$ and univariate normal base distribution,
$G_{0}=N\left(\mu_{G},\tau_{G}^{2}\right)$. The stick-breaking representation
of the Dirichlet process implies that distribution $G$ is almost surely
discrete because it has an infinite mixture distribution:

$G\stackrel{{\scriptstyle
d}}{{=}}\sum_{v=1}^{\infty}\varpi_{v}\delta_{\zeta_{v}},\text{ where
$\sum_{v=1}^{\infty}\varpi_{v}=1$ and $\zeta_{v}\stackrel{{\scriptstyle
iid}}{{\sim}}G_{0}$}.$ (5)

The distribution of random probabilities, $\varpi_{v}$, which depend on mass
parameter $\beta$, was derived in Sethuraman (1994); see also Ishwaran & James
(2003) and Lijoi & Prünster (2010). In the sequel, we condition on
distribution $G$; equivalently on the probabilities, $\varpi_{v}$, and
univariate atoms, $\zeta_{v}$, for $v\in\mathcal{N}$.

##### Cuisine 1 menu.

Recall that cuisine 1 represents the non-differential state, for which the
$T$-variate random vectors (i.e., dishes in the metaphor) has all equal
elements. Cuisine 1 menu, with its countably infinite dishes and their
associated probabilities, is modeled as follows as a discrete $T$-variate menu
distribution, $W_{1}$. With $\boldsymbol{1}_{T}$ denoting the column vector of
$T$ ones, let

$\boldsymbol{\vartheta}\mid\psi=\psi\boldsymbol{1}_{T}\text{ where $\psi\sim
G$.}$ (6)

Cuisine 1 menu distribution, $W_{1}$, is defined as the law of random vector
$\boldsymbol{\vartheta}$. Then $\mathcal{S}_{1}=$
$\\{\zeta_{v}\boldsymbol{1}_{T}:v\in\mathcal{N}\\}$ is the list of available
cuisine 1 dishes and the support of $W_{1}$. The selection probability
associated with dish $\zeta_{v}\boldsymbol{1}_{T}$ is $\varpi_{v}$.

The continuity of base distribution $G_{0}$ in equation (5) guarantees that
the menu $W_{1}$ dishes are unique. The discreteness of distribution $G$ has
practical implications for 2R2CF: (a) cuisine 1 consists of discrete dishes,
as required, rather than a continuous spectrum, and (b) since section 1 at
both the restaurants serve the same menu, the section 1 customers may eat the
same dish even if they select different restaurants.

##### Cuisine 2 menu.

As mentioned earlier, cuisine 2 depicts the differential state and its dishes
represent $T$-variate random vectors with at least two unequal elements. Its
menu comprises countably infinite cuisine 2 dishes along with their associated
probabilities. The menu is therefore modeled by a $T$-variate distribution,
$W_{2}$, satisfying two conditions: (i) it has a countably infinite support,
and (ii) each $T$-variate atom of $W_{2}$ has at least two unequal elements.
For any given
$\boldsymbol{\phi}=(\phi_{1},\ldots,\phi_{T})^{\prime}\in\mathcal{R}^{T}$, a
probability mass function for $W_{2}$ that satisfies these conditions is

$\displaystyle
W_{2}(\boldsymbol{\phi})=\begin{cases}\prod_{t=1}^{T}G(\phi_{t})/\bigl{(}1-\sum_{v=1}^{\infty}\varpi_{v}^{T}\bigr{)}&\qquad\text{if
$\phi_{t}\neq\phi_{t^{\prime}}$ for some $t\neq t^{\prime}$,}\\\
0&\qquad\text{otherwise,}\end{cases}$ (7)

where $G(\phi)$ denotes the mass function of distribution $G$, represented in
equation (5), evaluated at $\phi\in\mathcal{R}$. In line 1 of expression (7),
the normalizing constant $(1-\sum_{v=1}^{\infty}\varpi_{v}^{T})$ is the total
probability that a $T$-variate random vector whose elements are i.i.d. $G$ has
at least two distinct elements. Then menu distribution $W_{2}$ satisfies
conditions (i) and (ii), as required. Referring back to the atoms,
$\zeta_{v}$, of distribution $G$ in (5), the list of cuisine 2 dishes is
$\mathcal{S}_{2}=\\{\boldsymbol{\phi}:\phi_{t}=\zeta_{v}\text{ for some
}v\in\mathcal{N},\text{ but not all $\phi_{t}$ are identical},t=1,\ldots,T\\}$
and is the support of $W_{2}$. The selection probability associated with a
dish $\boldsymbol{\phi}$ is $W_{2}(\boldsymbol{\phi})$.

####  2.1.1 Restaurant, section, table, and dish choices of the 2R2CF
customers

Let the restaurant chosen by franchise customer $j$ be denoted by $g_{j}$ and
the chosen cuisine (i.e., section) be denoted by $s_{j}$. Suppose he or she
sits at table $v_{j}$ in that restaurant section and eats dish
$\boldsymbol{\theta}_{j}$.

##### Customer 1.

At time $0$, suppose the first customer selects restaurant $g_{1}=1$ with
probability $\rho_{1}>0$ and restaurant $g_{1}=2$ with positive probability
$\rho_{2}=1-\rho_{1}$. For reasons that will become clear, we refer to
$\rho_{1}$ as the baseline non-differential proportion and $\rho_{2}$ as the
baseline differential proportion. Typically, the differential state is
relatively less frequent, and so $\rho_{2}<\rho_{1}$ (i.e., $\rho_{2}<1/2$).
Proportion $\rho_{1}$ is given a uniform prior on the interval $(0.5,1)$.

Choice of cuisine $s_{1}$ Next, customer 1 selects a section in restaurant
$g_{1}$. As described earlier, each restaurant specializes in its namesake
cuisine, which is therefore more popular with its customers. This is modeled
as follows. Within restaurant $g_{j}$ (where $j=1$ for the first customer),
customer $j$ selects cuisine 1 with probability

$\displaystyle\mathcal{Q}_{g_{j}}(1)=\begin{cases}\rho_{1}+\rho_{2}\gamma&\qquad\text{if
$g_{j}=1$,}\\\ \rho_{1}-\rho_{1}\gamma&\qquad\text{if $g_{j}=2$,}\\\
\end{cases}$ (8)

for a speciality cuisine popularity parameter, $\gamma\in(0,1)$, that
determines the degree to which a restaurant’s patrons favor its namesake
cuisine. For instance, if $\gamma$ is nearly equal to 1, then a restaurant 1
(2) customer almost always (never) chooses cuisine 1. At the other extreme, if
$\gamma$ is nearly equal to 0, then the customer always chooses cuisine 1 with
approximate probability $\rho_{1}$ irrespective of the restaurant. Parameter
$\gamma$ is assigned an independent uniform prior on the unit interval. The
probability that a restaurant $g_{j}$ customer chooses cuisine 2 is then
$\mathcal{Q}_{g_{j}}(2)=1-\mathcal{Q}_{g_{j}}(1)$.

Choice of table $v_{1}$ and dish $\boldsymbol{\theta}_{1}$ Within section
$s_{1}$ of restaurant $g_{1}$, since the table identifiers are arbitrary, we
assume without loss of generality that customer 1 sits at table $v_{1}=1$. As
the 2R2CF process evolves, the tables in a restaurant’s section are assigned
consecutive labels as newly arriving patrons occupy them. At table $v_{1}=1$,
customer $1$ randomly orders a cuisine $s_{1}$ dish from menu distribution
$W_{s_{1}}$. The dish he or she eats represents the random effect of the first
probe. That is, $\boldsymbol{\theta}_{1}\mid s_{1}\sim W_{s_{1}}$.

Figure 2: Cartoon representation of the two-restaurant two-cuisine franchise
for differential analysis, showing the progressive choice of restaurant,
cuisine section, and table by customer $j$, for $j>1$. The numbered circles
represent table numbers. See the text in Section 2.1 for a detailed
description of the process.

##### Customer $j$, for $j>1$.

The restaurant choice of a subsequent customer is influenced by the previous
customer’s cuisine and the waiting time between the customers. Suppose
customer $j$ arrives at the franchise after a time interval of $e_{j-1}$
following the $(j-1)$th customer. Without loss of generality,
$e_{1},\ldots,e_{p-1}$ are scaled so that their total equals 1. Since the
probes in differential analysis typically represent CpG sites on a chromosome
or gene, it has a scaled length of 1.

To model the dependencies between the franchise customers, we first define a
non-negative dependence parameter $\eta$ that transforms the waiting time
$e_{j-1}$ to an affinity measure between customer $(j-1)$ and customer $j$:

$\displaystyle r_{j}=\exp(-e_{j-1}/\eta),\qquad j>1.$ (9)

The affinity measure belongs to the unit interval when $\eta>0$. If $\eta=0$,
affinity $r_{j}$ is defined as $0$ irrespective of the waiting time. The
affinity influences the behavior of customer $j$ through equation (10) below.

Choice of restaurant $g_{j}$ The cuisine $s_{j-1}$ of the $(j-1)$th customer
influences the restaurant choice of the $j$th customer through affinity
$r_{j}$ and popularity parameter $\gamma$:

$g_{j}\,\big{|}\,s_{j-1},\rho_{1},\eta,\gamma\sim\mathcal{F}_{j}.$

Specifically, the probability that customer $j$ selects restaurant 1 is
assumed to be

$\displaystyle\mathcal{F}_{j}(1)\stackrel{{\scriptstyle
def}}{{=}}P\bigl{(}g_{j}=1\mid
s_{j-1},\rho_{1},\eta,\gamma\bigr{)}=\begin{cases}\rho_{1}+\rho_{2}r_{j}/\gamma&\qquad\text{if
$s_{j-1}=1$,}\\\ \rho_{1}-\rho_{1}r_{j}/\gamma&\qquad\text{if $s_{j-1}=2$,}\\\
\end{cases}$ (10)

and $\mathcal{F}_{j}(2)=1-\mathcal{F}_{j}(1)$. If dependence parameter
$\eta=0$, we find that the restaurant choices of the customers are
independent; specifically, $\mathcal{F}_{j}(1)=\rho_{1}$ irrespective of the
cuisine $s_{j-1}$. The idea is illustrated in the top panel of Figure 2, where
customer $j$ chooses restaurant 1 with probability $\mathcal{F}_{j}(1)$ and
restaurant 2 with probability $\mathcal{F}_{j}(2)$.

It can be verified that $\mathcal{F}_{j}$ is a probability mass function if
and only if $r_{j}/\gamma<1$. Since the scaled waiting times are bounded above
by 1, a globally sufficient condition is $\eta<-1/\log\gamma$. We, therefore,
assume a mixture prior for dependence parameter $\eta$:

$\displaystyle\eta\mid\gamma\sim\frac{1}{2}\delta_{0}+\frac{1}{2}\mathcal{H}\cdot\mathcal{I}(\eta<-1/\log\gamma),$
(11)

where the second mixture component involves a continuous distribution,
$\mathcal{H}$, restricted to the interval $[0,-1/\log\gamma)$, thereby
enforcing the globally sufficient condition. In our experience, posterior
inferences on $\eta$ are relatively robust to the continuous prior
$\mathcal{H}$ provided the prior is not highly concentrated on a small part of
interval $[0,-1/\log\gamma)$.

When $\eta=0$, we have a zero-order Sticky PYP; when $\eta>0$, we obtain a
first order Sticky PYP. Some interesting consequences of specification (10)
are:

  1. 1.

Zero-order Sticky PYP: When $\eta=0$, each customer independently chooses
restaurant 1 (or 2) with a baseline probability of $\rho_{1}$ (or $\rho_{2}$).
The $p$ customers act identically.

  2. 2.

First order Sticky PYP with $e_{j-1}/\eta$ large: At large relative distances,
customer $j$ acts approximately independently of the history. Somewhat
similarly to customer 1, customer $j$ chooses restaurant 1 (2) with a
probability approximately, but not exactly, equal to $\rho_{1}$ ($\rho_{2}$).

  3. 3.

First order Sticky PYP with $e_{j-1}/\eta$ small: In the limit as
$e_{j-1}/\eta\to 0$ (e.g., for a small inter-probe distance $e_{j-1}$), the
restaurant choice of customer $j$ follows a hidden Markov model.

Since it drives the dependence characteristics of DNA methylation data,
parameter $\eta$ is of interest. Prior specification (11) allows the data to
direct the model order through posterior probability,
$P[\eta=0\mid\boldsymbol{X}]$, and an MCMC probability is readily available.

Choice of cuisine $s_{j}$ Upon entering restaurant $g_{j}$, customer $j$
selects cuisine-section $s_{j}$ with distribution $\mathcal{Q}_{g_{j}}$,
defined previously in expression (8). For bookkeeping reasons, among franchise
customers $1,\ldots,j$, suppose $p_{gs}^{(j)}$ customers choose section $s$ in
restaurant $g$; that is, $p_{gs}^{(j)}=$
$\sum_{l=1}^{j}\mathcal{I}(g_{l}=g,s_{l}=s)$ for $g,s=1,2$.

For a graphical depiction of cuisine selection by the $j$th customer, see the
middle panel of Figure 2, where $g_{j}=1$. That is, customer $j$, having
already chosen restaurant $1$, now chooses a cuisine-section. Restaurant 2 has
been greyed out because it is no longer accessible to this customer. In the
lower panel of Figure 2, we find that the customer picked cuisine-section 1,
so that $s_{j}=1$.

Choice of table $v_{j}$ Applying the above notation, among the previous
customers $1,\ldots,(j-1)$, there are $p_{g_{j}s_{j}}^{(j-1)}$ customers in
the same restaurant and section as the $j$th customer. Suppose these customers
have occupied tables $1,\ldots,M_{g_{j}s_{j}}^{(j-1)}$, and that there are
$p_{g_{j}s_{j}k}^{(j-1)}$ customers seated at the $k$th table. Let
$\mathcal{M}_{j-1}=\bigl{\\{}p_{g_{j}s_{j}k}^{(j-1)}:k=1,\ldots,M_{g_{j}s_{j}}^{(j-1)}\bigr{\\}}$
denote the aggregated table occupancies.

Recall that the newly arrived $j$th customer may sit at any of the
$M_{g_{j}s_{j}}^{(j-1)}$ occupied tables or a new
$(M_{g_{j}s_{j}}^{(j-1)}+1)$th table. Two possibilities are illustrated in the
lower panel of Figure 2. For a PYP with mass parameter $\alpha_{s_{j}}$ and
cuisine-specific discount parameter $d_{s_{j}}\in[0,1)$, the predictive
distribution of table $v_{j}$ of customer $j$ is assumed to be related to the
table occupancies:

$\displaystyle
P\biggl{(}v_{j}=k\mid\mathcal{M}_{j-1}\biggr{)}\propto\begin{cases}p_{g_{j}s_{j}k}^{(j-1)}-d_{s_{j}}\quad&\text{if
$k=1,\ldots,M_{g_{j}s_{j}}^{(j-1)}$,}\\\
\alpha_{s_{j}}+M_{g_{j}s_{j}}^{(j-1)}d_{s_{j}}\quad&\text{if
$k=(M_{g_{j}s_{j}}^{(j-1)}+1)$,}\\\ \end{cases}$ (12)

where the second line corresponds to customer $j$ sitting at a new table, in
which case the new number of occupied tables is
$M_{g_{j}s_{j}}^{(j)}=M_{g_{j}s_{j}}^{(j-1)}+1$ and table index
$v_{j}=M_{g_{j}s_{j}}^{(j)}$. Otherwise, if customer $j$ sits at a previously
occupied table, then their table index $v_{j}\leq M_{g_{j}s_{j}}^{(j-1)}$ and
the number of occupied tables remains unchanged:
$M_{g_{j}s_{j}}^{(j)}=M_{g_{j}s_{j}}^{(j-1)}$.

The above predictive distribution implies that customer $j$ is more likely to
choose tables with several occupants, positively reinforcing that table’s
popularity for future customers. The number of occupied tables stochastically
increases with the PYP mass and discount parameters.

For section $s$, if the PYP discount parameter $d_{s}=0$, we obtain the well-
known Pòlya urn scheme for Dirichlet processes (Ferguson, 1973b). PYPs act as
effective dimension reduction devices because the random number of occupied
tables is much smaller than the number of customers. In general, as the number
of patrons in section $s$ of restaurant $g$ grows as more customers arrive at
the franchise, that is, as $p_{gs}^{(j)}\to\infty$, the number of occupied
tables, $M_{gs}^{(j)}$, is asymptotically equivalent to

$\displaystyle\begin{cases}\alpha_{s}\log p_{gs}^{(j)}\qquad&\text{if
$d_{s}=0$}\\\ T_{d_{s}\alpha_{s}}\,p_{gs}^{(j)}\qquad&\text{if $0<d_{s}<1$}\\\
\end{cases}$ (13)

for a positive random variable $T_{d_{s}\alpha_{s}}$ (Lijoi & Prünster, 2010).
Asymptotically, the order of the number of occupied tables increases with
discount parameter $d_{s}$.

Choice of dish $\boldsymbol{\theta}_{j}$ As we have discussed earlier, the
2R2CF process assumes that all customers seated at a table of section $s_{j}$
are served the same dish, chosen from the cuisine $s_{j}$ menu by the first
customer to sit at that table. Let $\boldsymbol{\phi}_{g_{j}s_{j}k}$ denote
the common dish eaten by the previous $(j-1)$ customers at the $k$th table,
$k=1,\ldots,M_{g_{j}s_{j}}^{(j-1)}$.

Following predictive distribution (12), if customer $j$ has chosen to sit at a
previously occupied table, then she or he is served the dish selected by the
first customer to sit at that table. Otherwise, if customer $j$ has chosen a
new table, then she or he randomly picks a dish from menu distribution
$W_{s_{j}}$. The dish that customer $j$ eats represents the probe-specific
random effect $\boldsymbol{\theta}_{j}$, and

$\boldsymbol{\theta}_{j}\begin{cases}=\boldsymbol{\phi}_{g_{j}s_{j}v_{j}}\qquad&\text{if
$v_{j}=1,\ldots,M_{g_{j}s_{j}}^{(j-1)}$,}\\\ \sim W_{s_{j}}\qquad&\text{if
$v_{j}=(M_{g_{j}s_{j}}^{(j-1)}+1)$.}\\\ \end{cases}$ (14)

In the latter case (line 2), the dish $\boldsymbol{\theta}_{j}$ randomly
selected by customer $j$ is registered as
$\boldsymbol{\phi}_{g_{j}s_{j}M_{g_{j}s_{j}}^{(j)}}$ and is served to all
future customers who sit at table $M_{g_{j}s_{j}}^{(j)}$. Assumptions (12) and
(14) imply that although the restaurants serve the same menus, the relative
popularity of each dish is restaurant-specific.

The aforementioned process continues for the remaining 2R2CF customers.
Expressions (8) and (10) guarantee that a cuisine is more popular at its
namesake restaurant and the cuisine selected by a customer influences the
restaurant choice of the next customer, making the next customer likely to
select the same cuisine. This accounts for the lengthy runs of differential or
non-differential probes observed in methylation data. In addition to achieving
dimension reduction, the proposed Sticky PYP models the serial dependencies of
adjacent probes as a decreasing function of the inter-probe distances.

####  2.1.2 Latent clusters and their differential states

Latent clusters, introduced earlier in expression (3), comprise probes with
identical random effects and form the basis of the dimension reduction
strategy. Returning to the 2R2CF metaphor, we identify a cluster as the set of
customers who eat the same dish. However, in addition to the customers seated
at a table, multiple tables in both restaurants may serve the same dish
because of the shared cuisine menu. Therefore, irrespective of the restaurant,
aggregating customers eating the same dishes, we obtain the probe-cluster
allocation variables $c_{1},\ldots,c_{p}$, and the number of latent clusters,
$q$. The sets of customers eating the same cuisine 2 (differential state)
dishes correspond to the differential clusters, $\mathcal{D}$, defined in
equation (4).

From expression (13), we expect the number of occupied tables to be much
smaller than the number of franchise customers, $p$. Furthermore, since
multiple tables may serve the same dish, we expect the number of latent
clusters, $q$, to be smaller than the number of occupied tables. With high
probability, this implies that $q$ is much smaller than $p$.

##### PYP discount parameter $d_{2}$.

Consider the differential state cuisine menu, $W_{2}$, defined in (7). It can
be shown that as the number of treatments, $T$, and the number of probes, $p$,
increase, the differential clusters are not only asymptotically identifiable
but consistently detectable in the posterior; refer to Section 4 of Guha &
Baladandayuthapani (2016) for a detailed general discussion of this remarkable
phenomenon in standard PYP settings. Since the differential clusters can be
inferred with high accuracy when $T$ and $p$ are large, discount parameter
$d_{2}$ is given the mixture prior:

$\displaystyle d_{2}\sim\frac{1}{2}\delta_{0}+\frac{1}{2}U(0,1)$

where $d_{2}=0$ corresponds to a Dirichlet process. This provides the
posterior flexibility to choose between a Dirichlet process and a more general
PYP for a suitable clustering pattern of the differential probes. Allocation
patterns typical of Dirichlet processes, such as exponentially decaying
cluster sizes dominated by a few large clusters, result in high posterior
probabilities that $d_{2}$ equals $0$. By contrast, allocation patterns
characteristic of non-Dirichlet PYPs, such as slower-than-exponential power
law decays in the cluster sizes and relatively large numbers of smaller-sized
clusters, cause the posterior of discount parameter $d_{2}$ to concentrate
near 1 and exclude 0. A proof of the intrinsically different cluster patterns
of Dirichlet processes and PYPs is given in Theorem 2.1 of Guha &
Baladandayuthapani (2016).

Since distribution $G$ is discrete, all atoms of $T$-variate distribution
$W_{2}$ may not be unique. Indeed, this is common for $T=2$ treatments.
However, as $T$ grows, and provided the number of probes, $p$, grows at a
slower-than-exponential rate as $T$, the probability that two atoms allocated
to the probes are identical rapidly decays to 0. In regression problems
unrelated to differential analysis, Section 2.3 of Guha & Baladandayuthapani
(2016) derived a similar result for a simpler zero-order stochastic process.
We have verified this phenomenon in simulation studies on differential
analysis datasets. In several hundred artificial datasets generated from the
Sticky PYP, for $p=1,500$ probes and $T$ as small as four, no two allocated
atoms of $W_{2}$ were identical.

##### PYP discount parameter $d_{1}$.

Consider the (non-differential) cuisine menu, $W_{1}$, defined in (6). The
flexibility provided by PYP allocation patterns is not necessary for the non-
differential probes. This is because the allocation patterns of distribution
$W_{1}$ are driven by univariate parameter $\psi$ in (6), and, in general, the
mixture allocations of univariate objects are unidentifiable (e.g., Frühwirth-
Schnatter, 2006). Consequently, we set PYP discount parameter $d_{1}=0$,
reducing the two PYPs associated with the non-differential state (i.e.,
section 1 in both restaurants) to Dirichlet processes.

####  2.1.3 Other model parameters

Depending on the specifics of the application, an appropriate model is assumed
for the subject-specific parameters $\xi_{1},\ldots,\xi_{n}$. For example, we
may assume that $\xi_{i}\stackrel{{\scriptstyle
iid}}{{\sim}}N(0,\tau_{\epsilon}^{2})$. In some applications, it is more
appropriate to assume non-zero means and flexible error distributions:
$\xi_{i}=b_{i}+\epsilon_{i}$, where $b_{i}$ represents lane or batch effects
in methylation data, and the i.i.d. $\epsilon_{i}$ have a random distribution
with a Dirichlet process prior whose normal base distribution has mean zero
and variance $\tau_{\epsilon}^{2}$. Similarly, appropriate models for the
probe-specific parameters $\chi_{1},\ldots,\chi_{p}$ may include i.i.d. zero-
mean normal distributions, and finite mixtures or HMMs with state-specific
normal distributions. Inverse-gamma priors are assigned to $\sigma^{2}$ and
$\tau_{\epsilon}^{2}$. Appropriate priors are assumed for mass parameters
$\beta$, $\alpha_{1}$, and $\alpha_{2}$ in expressions (6) and (12). Mean
$\mu_{G}$ and variance $\tau_{G}^{2}$ of base distribution $G_{0}$ in
expression (6) are given a joint normal-inverse gamma prior. Figure 11 of
Supplementary Materials displays a directed acyclic graph summarizing the
relationships between the different BayesDiff model parameters.

##  3 Posterior Inference

Due to the analytical intractability of the BayesDiff model, we rely on MCMC
methods for posterior inferences and detection of differential probes.

###  3.1 MCMC Strategy

The model parameters are initialized using naïve estimation techniques and
iteratively updated by MCMC techniques until the chain converges. We split the
MCMC updates into three blocks. An outline of the MCMC procedure is as
follows. Further details can be found in Section 7 of Supplementary Material.

  1. 1.

Restaurant-cuisine-table-dish choice
$(g_{j},s_{j},v_{j},\boldsymbol{\theta}_{j})$ of customer $j$: For each probe
$j=1,\ldots,p,$ we sample the 4-tuple
$(g_{j},s_{j},v_{j},\boldsymbol{\theta}_{j})$ given the 4-tuples of the other
$(p-1)$ probes. This is achieved by proposing a new value of
$(g_{j},s_{j},v_{j},\boldsymbol{\theta}_{j})$ from a carefully constructed
approximation to its full conditional, and by accepting or rejecting the
proposed parameters in a Metropolis-Hastings step. The procedure is repeated
for all $p$ probes. As discussed in Section 2.1.2, the probe-cluster
allocations $c_{1},\ldots,c_{p}$ are immediately available from the
restaurant-cuisine-table allocations $(g_{j},s_{j},v_{j})$ of the $p$ probes.
Also available are the $q$ latent clusters along with their allocated probes,
and the set of differential clusters $\mathcal{D}$.

  2. 2.

Latent vectors $\boldsymbol{\lambda}_{1},\dots,\boldsymbol{\lambda}_{q}$:
There are $Tq$ latent vector elements, not all of which are necessarily
distinct because of the Dirichlet process prior on distribution $G$. Although
the latent vector elements are known from the Block 1 updates, MCMC mixing is
considerably improved by updating the latent vector elements conditional on
the $p$ probe-cluster allocations. As the calculation in Supplementary
Material shows, this is possible by Gibbs sampling.

  3. 3.

Remaining model parameters: Generated by standard MCMC techniques.

For the numerical analyses of this paper, we discarded an initial burn-in of
10,000 MCMC samples and used the subsequent 50,000 draws for posterior
inferences. Convergence was informally assessed by trace plots of various
hyperparameters to validate the MCMC sample sizes. For the proposed moves (in
discrete parameter space) described in Step 1, the average Metropolis-Hastings
acceptance rate exceeded 90% in all our analyses.

###  3.2 Detection of Differential Probes with FDR Control

Post-processing the MCMC sample, a Bayesian approach for controlling the false
discovery rate (FDR) (Newton et al., 2004) is applied to accurately detect the
probes $j$ with differential state $s_{j}=2$. Specifically, let $q_{0}$ be the
nominal FDR level and $\omega_{j}$ be the posterior probability that probe $j$
is differential, so that $\omega_{j}=P[s_{j}=2\mid\boldsymbol{X}]$. An
empirical average estimate, $\hat{\omega}_{j}$, is available from the MCMC
sample. To achieve the desired FDR level in calling the differential probes,
we first rank all the probes in decreasing order of $\hat{\omega}_{j}$. Let
$\hat{\omega}_{(1)}>\hat{\omega}_{(2)}>\dots>\hat{\omega}_{(p)}$ denote the
ordered posterior probability estimates. For each $b=1,\dots,p$, we calculate
as follows the posterior expected FDR resulting from calling the first $b$
probes in decreasing order of $\hat{\omega}_{j}$:

$\displaystyle\widehat{\text{FDR}}_{b}=\frac{\sum_{j=1}^{p}\left(1-\hat{\omega}_{j}\right)\mathcal{I}\left(\hat{\omega}_{j}\geq\hat{\omega}_{(b)}\right)}{\sum_{j=1}^{p}\mathcal{I}\left(\hat{\omega}_{j}\geq\hat{\omega}_{(b)}\right)}=\frac{\sum_{j=1}^{b}\left(1-\hat{\omega}_{(j)}\right)}{b},$
(15)

where the simplification occurs because the $\hat{\omega}_{j}$’s are sorted.
Finally, we pick the largest value of $b$, denoted by $b^{*}$, for which
$\widehat{\text{FDR}}_{b}<q_{0}$. A nominal FDR level of $q_{0}$ is achieved
by labeling as differential the first $b^{*}$ probes arranged in decreasing
order of $\hat{\omega}_{j}$.

##  4 Simulation Studies

Using artificial datasets with $T=5$ treatments, we analyzed the accuracy of
BayesDiff in detecting differentially methylated probes. We also compared the
results with established differential methylation procedures and general
statistical techniques for multiple treatment comparisons. Further, we
evaluated the ability of the BayesDiff procedure in discovering the complex
dependence structure of DNA methylation data.

##### Generation strategy

Proportions representing DNA methylation data were generated using the logit
transformation as in equation (1). The inter-probe distances were the actual
distances from the motivating TCGA dataset, scaled to add to $1$. In order to
capture the complexity of methylation data, such as the existence of multiple
latent methylation states (e.g,. CpG islands and shores), different read
depths across CpGs, and the incomplete conversion of bisulphite sequencing,
the generation strategy was partly based on techniques implemented in
WGBSSuite, a flexible stochastic simulation tool for generating single-base
resolution methylation data (Rackham et al., 2015). However, the generation
procedure differed from WGBSSuite in some respects. Specifically, it allowed
more than two treatments ($T>2$). Additionally, as in actual methylation
datasets, the generation procedure incorporated serial dependence not only in
the methylation levels but also in the differential states of the probes.

The probes-specific read depths were generated as
$n_{j}\stackrel{{\scriptstyle iid}}{{\sim}}\text{Poisson}(50)$. Unlike
BayesDiff assumption (1), there were no subject-specific random effects in the
generation mechanism. Instead, the normal mean of the generated data included
additive probe-specific random effects,
$\chi_{1}^{(0)},\ldots,\chi_{p}^{(0)}$, that were generated as follows:

  1. 1.

Generate the true methylation status of the probes, denoted by
$h_{1}^{(0)},\ldots,h_{p}^{(0)}$, using the 4-state distance-based HMM of
Rackham et al. (2015), with the states respectively representing the
methylated, first transit, demethylated, and second transit states.

  2. 2.

Set the baseline methylation levels for the methylated, (first or second)
transit, and demethylated states as $p_{\text{methylated}}=0.8$,
$p_{\text{transit}}=0.5$, and $p_{\text{un-methylated}}=0.2$.

  3. 3.

For $j=1,\ldots,p$, define the mean probe-specific random effect as follows:

$\displaystyle\tilde{\chi}_{j}^{(0)}=\begin{cases}\log(\frac{p_{\text{methylated}}}{(1-p_{\text{methylated}}})&\text{if
$h_{j}=1$ (i.e.,~{}methylated state),}\\\
\log(\frac{p_{\text{transit}}}{1-p_{\text{transit}}})&\text{if $h_{j}=2,4$
(first or second transit state),}\\\
\log(\frac{p_{\text{demethylated}}}{1-p_{\text{demethylated}}})&\text{if
$h_{j}=3$ (i.e.,~{}demethylated state).}\\\ \end{cases}$

  4. 4.

Independently generate true probe-specific random effects: $\chi_{j}^{(0)}\sim
N\bigl{(}\tilde{\chi}_{j}^{(0)},\tau_{\chi}^{2}\bigr{)}$ for probe
$j=1,\ldots,p$.

##### Noise and dependence levels

We investigated four scenarios corresponding to the combinations of two noise
levels and two dependence levels. For each scenario, 20 datasets were
independently generated, with each dataset consisting of $p=500$ probes and
$T=5$ treatments with 4 samples each, i.e., a total to $n=20$ samples. The low
noise level corresponded to true variance parameter $\sigma^{2}_{0}=0.36$;
equivalently, to a signal-to-noise of $R_{0}^{2}\approx 70\%$. The high noise
level corresponded to $\sigma^{2}_{0}=1$ or $R_{0}^{2}\approx 40\%$. The true
between-probe dependencies comprised two levels: no serial correlation (i.e.,
a zero-order Sticky PYP) with $\eta_{0}=0$, and positive serial correlation
(i.e., a first order Sticky PYP) with $\eta_{0}=0.004$. Although
$\eta_{0}=0.004$ may appear to be small, its value is calibrated to the inter-
probe distances; when the distance between two adjacent probes is equal to the
standardized average distance of $\bar{e}=1/(p-1)=1/499$, $\eta_{0}=0.004$
gives an affiliation of $r_{0}=0.6$ in equation (9). Since the affiliations
are bounded above by 1, $\eta_{0}=0.004$ represents fairly high inter-probe
dependence. For convenience, we will refer to the two dependence levels as
“no-correlation” and “high correlation.” The other model parameters were
common for the four scenarios and are displayed in Table 1. Setting a true
baseline differential proportion of $\rho_{2}=0.1$ resulted in approximately
10% true differentially methylated CpGs in each dataset.

$\alpha_{1}$ | $\alpha_{2}$ | $d_{2}$ | $\beta$ | $\gamma$ | $\rho_{2}$ | $\mu_{G}$ | $\tau^{2}_{G}$ | $\tau^{2}_{\chi}$
---|---|---|---|---|---|---|---|---
20 | 20 | 0.33 | 20 | 0.9 | 0.1 | 0 | 1 | 0.1225
Table 1: True parameter values used to generate the artificial datasets.

##### Posterior inferences

Assuming all model parameters to be unknown, each artificial dataset was
analyzed using a BayesDiff model that differed in key respects from the true
generation mechanism. For example, unlike the 4-state HMM of the generation
strategy, the probe-specific random effects $\chi_{j}$ were analyzed using a
BayesDiff model that ignored the first order dependence, and instead, relied
on a 3-state finite mixture model representing the methylated, transit, and
unmethylated states. Additionally, in contrast to the zeroed-out subject-
specific random effects during data generation, the BayesDiff procedure
assumed that the random effects were i.i.d. normal with zero mean.

To assess the accuracy of BayesDiff in detecting the absence or presence of
inter-probe serial correlation, in the no-correlation ($\eta_{0}=0$)
situation, we evaluated
$\log\left(\frac{P\left[\eta=0\mid\mathbf{X}\right]}{P\left[\eta>0\mid\mathbf{X}\right]}\right)$,
the log-Bayes factor comparing zero order to first order Sticky PYPs. In the
high correlation ($\eta_{0}=004$) situation, we evaluated
$\log\left(\frac{P\left[\eta>0\mid\mathbf{X}\right]}{P\left[\eta=0\mid\mathbf{X}\right]}\right)$,
the log-Bayes factor comparing first order to zero order Sticky PYPs. Thus, in
any scenario, a large positive value of this measure constitutes strong
evidence that BayesDiff detects the correct model order.

Although conceptually straightforward, the estimation of Bayes factors
requires multiple MCMC runs even for relatively simple parametric models
(Chib, 1995). Basu & Chib (2003) extended the estimation strategy to infinite
dimensional models such as Dirichlet processes. However, the computational
costs are prohibitively high for big datasets, and multiple MCMC runs for
estimating Bayes factors would stretch available computational resources
beyond their present-day limits. Faced with these challenges, we relied on an
alternative strategy for estimating the lower bounds of log-Bayes factors
using a single MCMC run. As it turns out, this is often sufficient to infer
the Sticky PYP model orders. Let $\Theta^{-}$ denote all BayesDiff model
parameters except $\eta$. In the high correlation situation, applying Jensen’s
inequality, a lower bound for the corresponding log-Bayes factor is
$E\left[\log\left(\frac{P\left[\eta>0\mid\mathbf{X},\Theta^{-}\right]}{P\left[\eta=0\mid\mathbf{X},\Theta^{-}\right]}\right)\mid\mathbf{X}\right]$.
Unlike log-Bayes factors, this lower bound can be easily estimated by an
empirical average estimate based on a single MCMC run. In the no-correlation
situation, a lower bound for the log-Bayes factor,
$\log\left(\frac{P\left[\eta=0\mid\mathbf{X}\right]}{P\left[\eta>0\mid\mathbf{X}\right]}\right)$,
can be similarly derived.

Figure 3: Simulation study box plots for estimated lower bounds of log-Bayes
factors in favor of the true model order.

For each of the four generation scenarios, box plots of these estimated lower
bounds for the 20 datasets are depicted in Figure 3. Except for the high
noise–no-correlation scenario, for which the results were inconclusive, the
estimated lower bounds of the log-Bayes factors in favor of the true
correlation structure were all positive and large. In the low noise–no-
correlation scenario, BayesDiff decisively favored zero-order models, and the
smallest lower bound among the 20 datasets was $13.9$, corresponding to Bayes
factors exceeding $e^{13.9}=1,088,161$. The $25^{th}$ percentile of these
lower bounds was 30.9, corresponding to Bayes factors exceeding
$e^{30.9}=2.63\times 10^{13}$. This is strong evidence that the BayesDiff
approach is reliable in this scenario. For the high-correlation scenarios, the
estimated lower bounds were even higher, showing that BayesDiff overwhelmingly
favors first order models when the data are actually serially correlated.

##### Comparisons with other methods

We evaluated the success of the BayesDiff procedure in detecting disease
genomic signatures and made comparisons with six well-known procedures. These
included some general statistical techniques for multiple treatment
comparisons, namely, one-way analysis of variance (ANOVA) and the Kruskal-
Wallis test. We also made comparisons with some methods specially developed
for detecting differential methylation in more than two treatments: COHCAP
(Warden et al., 2013), methylKit (Akalin et al., 2012), BiSeq (Hebestreit et
al., 2013), and RADMeth (Dolzhenko & Smith, 2014). The ANOVA and Kruskal-
Wallis test procedures were applied separately on each probe after applying
the inverse-logit transform. Being specifically designed for differential
methylation analysis, the COHCAP method was directly applied to the generated
proportions of the synthetic data. The remaining three methods are designed
for bisulfite sequencing, which consists of the total methylation reads for
each measured CpG site. For these methods, the methylation reads for each
probe was obtained by multiplying the proportion methylation values by the
total read. The bandwidth smoothing parameter of the method BiSeq was tuned to
optimize the overall detection. For all six methods, probe-specific p-values
were evaluated and adjusted for multiplicity using the FDR control procedure
of Benjamini & Hochberg (1995).

As with most nonparametric Bayes techniques, the associated computational
costs of BayesDiff are considerably higher compared to frequentist methods.
However, the additional run times are negligible compared to the time frames
over which the experimental data are collected, and as we demonstrate, the
trade-off is the substantially greater accuracy achieved by BayesDiff. On a
personal computer with an Intel Core i7-4770 processor with 3.40 GHz frequency
and 8 GB RAM, the average computational time for the Section 3.1 MCMC
algorithm, for the $n=20$ samples, $T=5$ treatments, and $p=500$ probes of the
simulation study datasets, was 0.60 seconds per iteration. The computational
times were reduced when different simulation datasets were run in parallel
across multiple cores of a research computing cluster. Analyzing datasets of
different sizes, we find that the computational cost is $O(Tp^{2})$ but does
not appreciably depend on $n/T$. This is reasonable because the mixture model
primarily focuses on
$\boldsymbol{\theta}_{1},\ldots,\boldsymbol{\theta}_{p}\in\mathcal{R}^{T}$.
Due to the intensive nature of the one-parameter-at-a-time Gibbs sampling
updates in Block 2, the Metropolis-Hastings algorithm of Guha (2010) can be
applied to significantly speed up the updates and make the calculations more
scalable. As part of ongoing development of a fast R package, we find that
ten- to hundred-fold speedups are possible with this fast MCMC strategy, which
can also be applied to quickly generate the Block 1 parameters.

We computed the receiver operating characteristic (ROC) curves for
differential probe detection for all seven methods. For a quantitative
assessment, we calculated the area under curve (AUC), declaring the method
with the largest AUC as the most reliable in each scenario. The ROC curves,
averaged over the 20 datasets under each simulation scenario, are shown with
the AUCs in Figure 7 of Supplementary Material. In all except the high-
noise–no-correlation scenario, BayesDiff uniformly outperformed the other
methods. Even in the high-noise–no-correlation scenario, BayesDiff performed
better in the low FPR region. As expected, all seven methods had lower
accuracies for higher noise levels. BayesDiff did significantly better than
the competing methods in the high correlation scenarios, suggesting that the
incorporation of between-probe dependencies improves its accuracy in
situations typical of DNA methylation data.

Since researchers typically focus on small false positive rates (FPRs), that
is, small significance levels, we also calculated the measures,
$\textrm{AUC}_{20}$ and $\textrm{AUC}_{10}$. $\textrm{AUC}_{20}$
($\textrm{AUC}_{10}$) is defined as the area under the ROC curve multiplied by
$5$ ($10$) when the FPR does not exceed $0.2$ ($0.1$). The multiplicative
factor ensures that the areas potentially vary between $0$ and $1$. The three
versions of AUC are presented in Table 3 in Supplementary Material. As also
seen in Figure 7, Table 3 reveals that in three of the four scenarios,
BayesDiff had the largest overall AUC. Furthermore, BayesDiff had vastly
improved reliability for low FPRs. For example, consider the low noise–high
correlation scenario. The overall AUC for BayesDiff was 0.035 greater than
that for ANOVA. In contrast, the gains for BayesDiff, relative to ANOVA, were
$+0.107$ for $\textrm{AUC}_{20}$ and $+0.146$ for $\textrm{AUC}_{10}$. The
advantages of BayesDiff were even greater relative to the other competing
methods. In the high noise–low-correlation scenario, BayesDiff had a
relatively low AUC, as mentioned. However, even in this scenario, it had the
greatest $\textrm{AUC}_{20}$ and $\textrm{AUC}_{10}$ values among all the
methods. Additionally, for a nominal FDR level of $q_{0}=0.05$, the achieved
FDR of BayesDiff was between 0 and 0.03 in every dataset and simulation
scenario. These results demonstrate the ability of BayesDiff to accurately
detect the differential probes even in challenging situations where the FPR is
small.

##  5 Data Analysis

We returned to the motivating DNA methylation dataset consisting of the upper
GI cancers: stomach adenocarcinoma (STAD), liver hepatocellular carcinoma
(LIHC), esophageal carcinoma (ESCA) and pancreatic adenocarcinoma (PAAD).
Applying the BayesDiff procedure, we detected the differentially methylated
CpG loci among the cancer types.

##### Data processing

The dataset was obtained from The Cancer Genome Atlas project, publicly
available through The Genomic Data Commons (GDC) portal (Grossman et al.,
2016). The data are available from the Illumina Human Methylation 450 platform
for each of 485,577 probes at the CpG sites. At the time of analysis, the
dataset consisted of $1,224$ tumor samples. The analysis was performed on a
gene-by-gene basis. We picked a set of 443 genes involving mutation in at
least 5% of the samples. To ensure that all CpG sites potentially related to a
gene were included in the analysis, we selected all sites located within 50K
base pairs outside the gene body, specifically, upstream from the 5’ end and
downstream from the 3’ end. The number of gene-specific CpG sites ranged from
1 to 769, and are displayed in Figure 8(a) of Supplementary Material. As a
final preprocessing step, since the methylation patterns of short genes are of
less interest in cancer investigations, we eliminated the 25 genes mapped to
20 or fewer CpG sites.

##### Inference procedure

The data were analyzed using the proposed BayesDiff approach. Exploratory
analyses indicated that a satisfactory fit was obtained by eliminating the
probes-specific random effects $\chi_{j}$ in (1). Since information about
experimental batches is not available for the TCGA dataset, we assumed that
the i.i.d. $\xi_{1},\ldots,\xi_{n}$ parameters follow a random distribution
with a zero-mean Dirichlet process prior. The MCMC procedure of Section 3.1
was applied to obtain posterior samples for each gene. For detecting
differential CpG sites, we applied the Section 3.2 procedure with a nominal
FDR of $q_{0}=0.05$.

##### Results

(a) Contingency table of detected methylation status and location of CpG site
with respect to gene body

(b) Contingency table of detected methylation status and proximity of CpG site
to chromosomal end

Figure 4: Associations of detected methylation status and position of CpG
sites.

Among the differentially methylated CpG sites detected by our approach,
approximately 40.6% of the sites were located outside the gene bodies. Figure
4 displays the associations between detected methylation status and positions
of the CpG sites. For our analysis, we have defined “near the 5’ (3’) end” as
the CpG sites located within one-fourth length of the gene body, either inside
or outside the gene boundary, and closer to the transcription start
(termination) site. Our results indicate that the proportion of differential
methylation is higher for CpG sites inside the gene body, and that most
differentially methylated loci are situated within the gene body, as is well
known from numerous previous studies. However, our analysis also revealed
significant amounts of differential methylation outside the gene body. Despite
the common belief that DNA methylation analysis should focus on the 5’ end
region, we found that CpG sites near the 3’ ends also displayed considerable
degrees of differential methylation. These findings support the
recommendations of Irizarry et al. (2009) that studies of DNA methylation
alteration should be conducted on a higher resolution, epigenome-wide basis.

Figure 5: Site-wise summary of the largest pairwise differences of
differentially methylated loci among the four upper GI cancer types.

Among the differentially methylated sites detected by BayesDiff, we estimated
the pairwise differences between random effects associated with the four
cancer types. Site-wise summaries of the largest pairwise differences of the
cancer-specific effects are displayed in Figure 5. None of the four cancer
types displayed consistent hypermethylation or hypomethylation across all
genes or entire chromosomes. However, we found that LIHC is frequently
differentially methylated relative to one of the other cancer types, implying
that it is the most volatile disease with respect to DNA methylation.

For each gene, Figure 8(b) of Supplementary Material displays 95% credible
intervals for lower bounds of log-Bayes factors of a first versus zero-order
model, i.e., $\eta=0$ versus $\eta>0$ in expression (9). Models with first
order dependence are overwhelmingly favored for a majority of the genes,
suggesting that statistical techniques that do not account for dependence
between neighboring CpG sites would be less effective for these data. Figure 9
of Supplementary Material displays the detailed differential methylation
pattern for the top two mutated genes, TP53 and TTN. An obvious feature of
both the genes is that the detected differential methylation of the CpG sites
is highly serially correlated. For gene TP53, there are almost no
differentially methylated loci within the gene body. The 3’ end region outside
the gene body has a cluster of differentially methylated loci, for which
cancer type STAD is mostly hypermethylated. The results for gene TTN tell a
quite different story: most of the differentially methylated loci are inside
the gene body and near the 5’ end. Cancer type LIHC is hypomethylated compared
to PAAD around the 5’ end region, but it is hypermethylated compared to STAD
near the 3’ end. Genes with at least 90% differentially methylated sites
detected by BayesDiff are listed in Table 4 of Supplementary Material, along
with the largest pairwise difference between the four cancer types among the
differentially methylated loci. The number of CpG sites within each segment is
listed in Table 5 of Supplementary Material.

Existing medical literature both supports and complements our findings. For
example, hypermethylation of the EDNRB and SLIT2 genes have been found in STAD
(Tao et al., 2012). Gene FBN2 was hypermethylated in ESCA (Tsunoda et al.,
2009). While several studies have found that the gene and protein expressions
of ABC transporter genes, such as ABCC9, are useful for understanding the
prognosis of esophageal cancer (Vrana et al., 2018), we conclude that
hypermethylation of ABCC9 is a major difference between cancer types ESCA and
LIHC. Gene FLRT2 is a potential tumor suppressor that is hypermethylated and
downregulated in breast cancer (Bae et al., 2017). Our results indicate that
this gene is also hypermethylated in cancer type STAD versus LIHC. Mutations
in SPTA1 gene has been linked with PAAD (Murphy et al., 2013). Our results
indicate that hypermethylation of this gene distinguishes PAAD from LIHC.

Finally, we compared our findings with those of ANOVA for multiple treatment
comparisons. Table 6 lists the common set of genes with at least 90%
differentially methylated sites identified by both BayesDiff and ANOVA. Table
7 displays the genes identified by only ANOVA, whereas Table 8 displays the
large number of genes detected exclusively by BayesDiff. Cross-referencing
with the medical literature, we find that genes FLRT2 and FBN2 were detected
by both methods. However, genes EDNRB, SLIT2, ABCC9, and SPTA1 were only
identified by BayesDiff, revealing the benefits of the proposed method.

##### Accounting for data characteristics

A statistical model for biomarker data should account for the observed probe-
specific means and variances. This is especially important in multiple-testing
based approaches where the first two sample moments must be plausibly
explained by the fitted model to avoid making misleading biological
interpretations (Subramaniam & Hsiao, 2012). From this perspective, certain
aspects of the BayesDiff model, such as variance $\sigma^{2}$ being a priori
unrelated to the mean in expression (1), may appear to be restrictive.
However, even though the BayesDiff model was not specifically constructed to
match data summaries such as sample moments, in practice, the nonparametric
nature of the Sticky PYP allows the posterior to flexibly adapt to the
features of the data, including sample moments, thereby accounting for mean-
variance relationships in a robust manner. For example, consider again the top
mutated genes, TP53 and TTN, discussed in Figure 9 of Supplementary Material.
The ability of the BayesDiff model to match the sample moments of the gene-
specific probes is demonstrated as follows. Given the inter-probe distances,
the joint posterior of the BayesDiff parameters induces predictive
distributions on the $n$ measurements for each probe. Functionals of these
predictive distributions, such as the probe-specific sample moments, are
easily estimated by post-processing the MCMC sample. For these two genes,
Figure 10 of Supplementary Material shows that the sample moments predicted by
BayesDiff closely match the actual first and second sample moments, with
correlations exceeding 99% in each plot. Similar results were observed in
other datasets.

##  6 Discussion

DNA methylation data exhibit complex structures due to unknown biological
mechanisms and distance-dependent serial correlations among neighboring CpG
sites or probes. The identification of differential signatures among multiple
treatments or sets of tumor samples is crucial for developing targeted
treatments for disease. This paper formulates a flexible approach applicable
to multiple treatments called BayesDiff. The technique relies on a novel
Bayesian mixture model called the Sticky PYP or the two-restaurant two-cuisine
franchise. In addition to allowing simultaneous inferences on the probes, the
model accommodates distance-based serial dependence and accounts for the
complex interaction patterns commonly observed in cancer data. A effective
MCMC strategy for detecting the differential probes is developed. The success
of the BayesDiff procedure in differential DNA methylation, relative to well-
established methodologies, is exhibited via simulation studies. The new
technique is applied to the motivating TCGA dataset to detect the differential
genomic signatures of four upper GI cancers. The results both support and
complement various known facts about the epigenomic differences between these
cancer types, while revealing a set of genes exhibiting high proportions of
differentially methylated CpG sites.

In addition to providing a good fit for the data, a statistical model must be
able to account for features of the dataset such as sample moments. The
success of the BayesDiff model in this regard is demonstrated in Section 5 for
the upper GI dataset. It must be emphasized, however, that there may be
datasets where BayesDiff is less successful in accounting for the data
characteristics, possibly due to slow asymptotic convergence of the posterior
to the true generative process. In such situations, more flexible global
transformations (Li et al., 2016) or variance-stabilizing transformations
(Durbin et al., 2002) may be used. Alternatively, local Laplace approximations
of exponential family likelihoods through link functions (Zeger & Karim, 1991;
Chib & Winkelmann, 2001) may extend the BayesDiff procedure and better explain
unique data characteristics.

Like most Bayesian models comprising several latent parameters, the proposed
2R2CF can be marginalized over different parameter sets to obtain equivalent
versions of the same model. For example, we could marginalize over the
restaurants to obtain an equivalent “sticky cuisine” version in which there is
just one restaurant with two cuisine-sections and a customer more likely to
favor the cuisine selected by the previous customer. Alternatively, we could
integrate out the sections to obtain an equivalent “sticky restaurant
franchise” in which each restaurant comprises a single section with
restaurant-specific probabilities that guarantee that Cuisine 1 or 2 dishes
are more popular at their namesake restaurant; a customer is then more likely
to favor the restaurant selected by the previous customer. In all equivalent
versions, however, a probe’s differential state is determined by the
characteristics of the customer’s dish in the metaphor.

The 2R2CF perspective offers the twin advantages of parameter interpretability
and generalizability. Section 9 of Supplementary Material presents the
generalized form of the Sticky PYP, revealing the full potential of the
proposed method in analyzing not only DNA methylation datasets, but also other
types of omics datasets, such as gene expression, RNASeq, and copy-number
alteration data. Beyond biomedical applications, the generalized formulation
offers a diverse palette of parametric and nonparametric models for capturing
the distinctive features of datasets. These Bayesian mixture models are
special cases of Sticky PYPs for particular choices of a countable group
parameter (e.g., two “restaurants” in the 2R2CF metaphor for differential
methylation problems) and countable state parameter (e.g., two “cuisines” in
2R2CF) with the state of a customer influencing the group of the next
customer. In addition to extending PYPs to discrete time series type data, the
range of models includes Dirichlet processes, PYPs, infinite HMMs,
hierarchical Dirichlet processes (Teh et al., 2006), hierarchical Pitman-Yor
processes (Teh et al., 2006; Camerlenghi et al., 2019), finite HMMs, nested
Chinese restaurant processes (Blei & Jordan, 2005), nested Dirichlet processes
(Rodriguez et al., 2008), and analysis of densities models (Tomlinson &
Escobar, 2003).

Ongoing work involves extending the correlation structure to model more
sophisticated forms of inter-probe dependence in DNA methylation data.
Commented R code implementing the BayesDiff method is available on GitHub at
https://github.com/cgz59/ BayesDiff. Using high-performance Rcpp subroutines,
we are developing a fast R package for detecting differential genomic
signatures in a wide variety of omics datasets. Initial results indicate that
dramatic speedups of several orders of magnitude would allow the fast analyses
of high-dimensional datasets on personal computers.

SUPPLEMENTARY MATERIAL

##  7 MCMC Strategy

As mentioned in Section 3.1 of the paper, we split the MCMC updates into three
blocks. The MCMC procedure for updating Blocks 1 and 2 is described below.
Unless otherwise stated, all references to equations, tables, and figures are
for the main paper.

  1. 1.

Restaurant-cuisine-table-dish choice
$(g_{j},s_{j},v_{j},\boldsymbol{\theta}_{j})$ of customer $j$: For each probe
$j=1,\ldots,p,$ we sample the 4-tuple
$(g_{j},s_{j},v_{j},\boldsymbol{\theta}_{j})$ given the 4-tuples of the other
$(p-1)$ probes. This is achieved by proposing a new value of
$(g_{j},s_{j},v_{j},\boldsymbol{\theta}_{j})$ from a carefully constructed
approximation to its full conditional, and by accepting or rejecting the
proposed 4-tuple of probe-specific parameters in a Metropolis-Hastings step.
The procedure is repeated for all $p$ probes to complete the block of MCMC
updates.

We give here an intuitive description of the details. Further calculation
details are available in Section 8 of Supplementary Material. For the $j$th
probe, $1<j<p$, the choice of restaurant $g_{j}$ depends on the triplets of
the two immediately adjacent probes. Specifically, as discussed in Section 2.1
and graphically shown in the upper and middle panel of Figure 2, the
restaurant selected by a customer depends on the cuisine of the previous
customer. Then, within restaurant $g_{j}$, the customer chooses cuisine
$s_{j}$ with a restaurant-specific probability. Finally, the customer may
either sit an existing table, in which case they must eat the same dish as
everyone already sitting at that table, or the customer may open a new table
and order a dish from cuisine menu $s_{j}$. The generation strategy for table
and dish depends on the cuisine, as discussed below.

##### Cuisine $s_{j}=1$

Evaluating the posterior probability of table $v_{j}$ involves integrating the
Gaussian likelihood of column vector
$\mathbf{z}_{j}=(z_{1j},\ldots,z_{nj})^{\prime}$ with respect to the cuisine 1
menu, i.e., $\boldsymbol{\theta}_{j}\sim W_{1}$. Recall that the cuisine 1
menu has the special structure (6), allowing its reduction to a univariate
quantity distributed as $G$. This random distribution itself follows a
Dirichlet process conditional prior with a base distribution that is a mixture
of known atoms (given the 4-tuples of the other probes) and Gaussian
distribution $G_{0}$. This conjugate hierarchical structure allows the
calculation of table and dish choice for cuisine 1 in computationally closed
form.

##### Cuisine $s_{j}=2$

This situation is more complicated because the Gaussian likelihood must be
integrated with respect to prior distribution $W_{2}$ of the cuisine 2 dishes.
Unlike cuisine 1, this is not possible in computationally closed form because
menu structure (7) cannot be reduced to a univariate quantity. In other words,
for cuisine 2 (i.e., the differential state), the Dirichlet process
conditional prior for distribution $G$ is non-conjugate.

We utilize an auxiliary variable approach. Given the current set and frequency
of atoms of distribution $G$ in the other $(p-1)$ probes, we first compute a
finite mixture approximation, $G^{*}$, to the Dirichlet process conditional
prior for $G$, as in Ishwaran & Zarepour (2002). Then, using prior (7), it is
possible to approximate the posterior probability of customer $j$ joining an
existing table or sitting at a new table, irrespective of the dishes. This
allows us to generate the table selection, $v_{j}$. Finally, we propose the
dish, $\boldsymbol{\theta}_{j}=\boldsymbol{\phi}_{g_{j}s_{j}v_{j}}$. If table
$v_{j}$ is one of the already existing tables, then the customer simply
chooses the common, table-specific dish. Otherwise, we sequentially generate
from the posterior the $T$ elements representing the new dish. This is done
using the corresponding treatment-specific elements of vector
$\mathbf{z}_{j}=(z_{1j},\ldots,z_{nj})^{\prime}$ and finite mixture $G^{*}$,
with the restriction that not all $T$ atoms of vector
$\boldsymbol{\theta}_{j}$ are equal. The Metropolis-Hasting step compensates
for any approximations and produces post–burn-in samples from the BayesDiff
model posterior.

##### Latent clusters

As discussed in Section 2.1.2, the probe-cluster allocations
$c_{1},\ldots,c_{p}$ are immediately available from the restaurant-cuisine-
table allocations $(g_{j},s_{j},v_{j})$ of the $p$ probes. The $q$ latent
clusters, along with their allocated probes, and the set of differential
clusters $\mathcal{D}$, are also immediately available.

  2. 2.

Latent vectors $\boldsymbol{\lambda}_{1},\dots,\boldsymbol{\lambda}_{q}$:
Since each latent vector is of length $T$, there are $Tq$ latent vector
elements, not all of which are necessarily distinct because of the Dirichlet
process prior on distribution $G$. Although the latent vector elements are
available from Block 1, mixing of the MCMC chain is considerably improved by
an additional block of updates of the latent vector elements conditional on
the $p$ probe-cluster allocations. As the following calculation shows, this
can be done by Gibbs sampling.

Since a non-differential cluster’s latent vector consists of a single value
repeated $T$ times, the number of distinct latent vector elements is no more
than $\left(q_{1}+Tq_{2}\right)$. Given the current probe-cluster allocations,
suppose there are $m_{k}$ probes associated with latent cluster $k$, where
$k=1,\dots,q$. Let $q_{2}=|\mathcal{D}|$ be the number of differential
clusters, so that $q_{1}=q-q_{2}$ is the number of non-differential clusters.

##### Non-differential clusters

For these clusters, the latent vector $\boldsymbol{\lambda}_{k}$ takes the
form, $\boldsymbol{\lambda}_{k}=\psi_{k}\boldsymbol{1}_{T}$ for some
$\psi_{k}\in\mathcal{R}$. Given the other $(q-1)$ latent vectors, assumptions
(6) and (7) imply that parameter $\psi_{k}$ has a Dirichlet process
conditional prior. The base distribution of the Dirichlet process is a mixture
of a continuous distribution, $G_{0}=N(\mu_{G},\tau_{G}^{2})$, and the
univariate atoms located at the “known” elements of the other $(q-1)$ latent
vectors. The sufficient statistic for $\psi_{k}$ is

$\displaystyle\bar{y}_{k}$
$\displaystyle=\frac{1}{nm_{k}}\sum_{j=1}^{p}\sum_{i=1}^{n}\left({z_{ij}}-{\xi_{i}}-\chi_{j}\right)\mathcal{I}\left(c_{j}=k\right)$
$\displaystyle\sim
N\left(\psi_{k},\frac{\sigma^{2}}{nm_{k}}\right)\quad\text{for cluster
}k\not\in\mathcal{D}.$

Since the Dirichlet process conditional prior is conjugate to this normal
likelihood, parameter $\psi_{k}$ could be updated by Gibbs sampling (Escobar,
1994; MacEachern, 1994; Escobar & West, 1995).

##### Differential clusters

For these clusters, some of the $T$ elements of latent vector
$\boldsymbol{\lambda}_{k}=(\lambda_{1k},\ldots,\lambda_{Tk})^{\prime}$ could
be tied, but at least two elements must be unequal. Denote these restrictions
on $\boldsymbol{\lambda}_{k}$ by $\mathcal{A}$. For a treatment $t$, let the
remaining $(T-1)$ elements of vector $\boldsymbol{\lambda}_{k}$ be denoted by
$\boldsymbol{\lambda}_{-tk}$.

Imagine updating element $\lambda_{tk}$ assuming that
$\boldsymbol{\lambda}_{-tk}$ and the remaining $(q-1)$ latent vectors are
known. Vector $\boldsymbol{\lambda}_{-tk}$ and restriction $\mathcal{A}$ imply
a possibly restricted support, $\mathcal{A}_{t}$, for parameter
$\lambda_{tk}$. Observe that the support is unrestricted, i.e.,
$\mathcal{A}_{t}=\mathcal{R}$, if at least two elements of
$\boldsymbol{\lambda}_{-tk}$ differ. Similarly to non-differential clusters,
it can be shown that element $\lambda_{tk}$ follows a Dirichlet process
conditional prior restricted to the support $\mathcal{A}_{t}$, and that the
base distribution is a mixture of a normal distribution and univariate atoms.
The sufficient statistic for latent vector element $\lambda_{tk}$ is

$\displaystyle\bar{y}_{tk}$
$\displaystyle=\frac{1}{n_{t}m_{k}}\sum_{j=1}^{p}\sum_{i:t_{i}=t}\left({z_{ij}}-{\xi_{i}}-\chi_{j}\right)\mathcal{I}\left(c_{j}=k\right)$
$\displaystyle\sim
N\left(\lambda_{tk},\frac{\sigma^{2}}{n_{t}m_{k}}\right)\quad\text{for
treatment $t=1,\ldots,T$, and cluster $k\in\mathcal{D}$},$

where $n_{t}$ is the number of individuals associated with treatment $t$. The
conjugate structure implies that parameter $\lambda_{tk}$ could be generated
from its full conditional by a rejection sampler on set $\mathcal{A}_{t}$
coupled with standard Gibbs proposals for conjugate Dirichlet processes. In
practice, the acceptance rates of the rejection sampler are very high for $T$
as small as 3 or 4, and are nearly 100% for larger $T$. This is because the
posterior probability of support set $\mathcal{A}_{t}$ tends to 1 as $T$
grows. An intuitive explanation for this asymptotic property is provided in
Section 2.1.3 following the prior specification of discount parameter $d_{2}$.

##  8 Additional details about Block 1 MCMC updates

In the first block of MCMC updates of Section 3.1, we sample the the set

$\mathcal{C}=\left\\{(g_{j},s_{j},v_{j},\boldsymbol{\theta}_{j}):j=1,\ldots,p\right\\}.$

We iteratively sample the 4-tuple,
$(g_{j},s_{j},v_{j},\boldsymbol{\theta}_{j})$, for each probe $j$ via
Metropolis-Hastings updates. Specifically, for the $j$th probe, we propose
group (i.e., restaurant) $g_{j}$, state (i.e., cuisine) $s_{j}$, PYP cluster
(i.e., table) label $v_{j}$, and random effect (i.e., dish)
$\boldsymbol{\theta}_{j}$, with the proposal distribution approximately equal
to the joint posterior of $(g_{j},s_{j},v_{j},\boldsymbol{\theta}_{j})$
conditional on the data and all other parameters being equal to the current
values. As explained in the paper, since the base measure of the Sticky PYP is
discrete and not conjugate to the likelihood, this makes the case of parameter
$v_{j}$ being assigned a new cluster label complicated. To deal with this
situation, an auxiliary variable approach is used. Conditional on the current
parameter values of the Dirichlet process prior for the latent vector
elements, we generate a finite-dimensional approximation, $G^{*}$, to the
Dirichlet process (Ishwaran & Zarepour, 2002). Since the distribution $G^{*}$
is discrete, it is characterized by a finite vector of probability masses,
$\boldsymbol{\pi}$, and the corresponding values of probability mass points,
$\mathbf{u}$. Then, conditional on the auxiliary variables from $G^{*}$, we
obtain the quantities needed for the conditional posterior of $v_{j}$ being
assigned a new cluster label in the Metropolis-Hastings algorithm.

More formally, consider probe $j\in\\{1,\dots,p\\}$. The proposal
probabilities for $\left(g_{j},s_{j},v_{j}\right)$ are as follows:

$\displaystyle\bullet\quad\text{For an existing PYP cluster indexed by
}v\in\\{1,\ldots,q_{gs}^{(-j)}\\},$ $\displaystyle
Q\left(g_{j}=g,s_{j}=s,v_{j}=v\right)=P\left(g_{j}=g,s_{j}=s,v_{j}=v\mid\mathbf{z}_{j},\mathbf{c}^{(-j)},\boldsymbol{\phi}_{gsv},s_{j-1},\eta,\rho,\gamma,\boldsymbol{\xi},\boldsymbol{\chi},\sigma^{2}\right)$
$\displaystyle\propto\left[\mathbf{z}_{j}\mid
v,s,g,\boldsymbol{\phi}_{gsv},\boldsymbol{\xi},\boldsymbol{\chi},\sigma^{2}\right]P\left(v_{j}=v\mid\mathbf{c}^{(-j)},s,g\right)P\left(s_{j}=s\mid
g,\gamma,\rho\right)P\left(g_{j}=g\mid s_{j-1},\eta,\gamma,\rho\right)$
$\displaystyle\propto\prod_{i=1}^{n}\varphi\left({z_{ij}}\mid\phi_{t_{i}gsv}+\xi_{i}+\chi_{j},\sigma^{2}\right)\left(m^{(-j)}_{gsv}-d_{s}\right)\mathcal{Q}_{g}(s)\mathcal{F}_{j}(g).$
(16) $\displaystyle\bullet\quad\text{For a new PYP cluster indexed by
}v^{*}=q_{gs}^{(-j)}+1,$ $\displaystyle
Q\left(g_{j}=g,s_{j}=s,v_{j}=v^{*}\right)=P\left(g_{j}=g,s_{j}=s,v_{j}=v^{*}\mid\mathbf{z}_{j},\mathbf{c}^{(-j)},\boldsymbol{\pi},\mathbf{u},s_{j-1},\eta,\rho,\gamma,\boldsymbol{\xi},\boldsymbol{\chi},\sigma^{2}\right)$
$\displaystyle\propto\left[\mathbf{z}_{j}\mid
v^{*},s,g,\boldsymbol{\pi},\mathbf{u},\boldsymbol{\xi},\boldsymbol{\chi},\sigma^{2}\right]P\left(v_{j}=v^{*}\mid\mathbf{c}^{(-j)},s,g\right)P\left(s_{j}=s\mid
g,\gamma,\rho\right)P\left(g_{j}=g\mid s_{j-1},\eta,\gamma,\rho\right)$
$\displaystyle\propto\begin{cases}\\!\begin{aligned}
&\sum_{l=1}^{L}\pi_{l}\left\\{\prod_{i=1}^{n}\varphi\left(z_{ij}\mid
u_{l}+\xi_{i}+\chi_{j},\sigma^{2}\right)\right\\}\left(\alpha_{1}+d_{1}q^{(-j)}_{g1}\right)\mathcal{Q}_{g}(1)\mathcal{F}_{j}(g),\quad\text{if
$s=1$,}\\\ \end{aligned}\\\ \\!\begin{aligned}
\left\\{\prod_{t=1}^{T}\sum_{l=1}^{L}\pi_{l}\prod_{i:t_{i}=t}\varphi\left({z_{ij}}\mid
u_{l}+{\xi_{i}}+\chi_{j},\sigma^{2}\right)-\sum_{l=1}^{L}\pi_{l}^{T}\prod_{i=1}^{n}\varphi\left({z_{ij}}\mid
u_{l}+{\xi_{i}},\sigma^{2}\right)\right\\}\\\
\times\left(\alpha_{2}+d_{2}q^{(-j)}_{g2}\right)\mathcal{Q}_{g}(2)\mathcal{F}_{j}(g),\quad\text{if
$s=2$.}\\\ \end{aligned}\\\ \end{cases}$ (17)

In expressions (16) and (17), when $j<p$, superscript $(-j)$ for a variable
indicates that the calculation excludes the $j$th and $(j+1)$th cluster
allocation variables. Specifically, $q^{(-j)}_{gs}$ denotes the number of PYP
clusters for group $g$ and state $s$; $\mathbf{c}^{(-j)}$ denotes the vector
of cluster allocations; $m^{(-j)}_{gsv}$ denotes the cluster membership count
for cluster $v$ in the PYP with group $g$ and state $s$; $\mathcal{Q}_{g}(s)$
and $\mathcal{F}_{j}(g)$ are defined in the paper;
$\varphi\left(\cdot\mid\mu,\sigma^{2}\right)$ is the density of the normal
distribution with mean $\mu$ and variance $\sigma^{2}$; $L$ is the number of
distinct mass points in $G^{*}$; $u_{l}$ denotes the value of distinct mass
points in $G^{*}$ and $\pi_{l}$ is the corresponding probability mass. Before
calculating the proposal, latent vectors of emptied clusters are removed and
cluster labels rearranged so that the largest label is $q^{(-j)}_{gs}$.

The case $v^{*}=(q_{gs}^{(-j)}+1)$ corresponds to the case of $v_{j}$ opening
a new cluster. If $v_{j}=(q_{gs}^{(-j)}+1)$ has been proposed, then a new
random effect of length $T$ is sampled from the posterior distribution based
on the $n$-dimensional observation vector $\mathbf{z}_{j}$ and the prior
$G^{*}$ characterized by the auxiliary variables, subject to the constraints
imposed on the random effects of the differential state $s_{j}$.

The Metropolis-Hastings acceptance ratio is computed to decide whether the
proposed $g_{j}$, $s_{j}$ and $v_{j}$ values are accepted. For $j=p$, our
proposal density is exactly the desired conditional posterior, so we always
accept the move. For $j<p$, since our proposal density is part of the desired
conditional posterior for $\left(g_{j},s_{j},v_{j}\right)$, this part cancels
out in the acceptance ratio. Therefore, the acceptance ratio only relates to
the transition from $\left(g_{j},s_{j},v_{j}\right)$ to
$\left(g_{j+1},s_{j+1},v_{j+1}\right)$, which is

$\displaystyle r_{j}=\frac{p\left(g_{j+1},s_{j+1},\tilde{v}_{j+1},\mid
g_{j}^{*},s_{j}^{*},v_{j}^{*},\boldsymbol{\phi}^{*}_{g_{j+1}s_{j+1}\tilde{v}_{j+1}},\mathbf{z}_{j+1},\mathbf{c}^{-},\eta,\boldsymbol{\xi},\boldsymbol{\chi},\sigma\right)}{p\left(g_{j+1},s_{j+1},v_{j+1}\mid
g_{j}^{0},s_{j}^{0},v_{j}^{0},\boldsymbol{\phi}_{g_{j+1}s_{j+1}v_{j+1}},\mathbf{z}_{j+1},\mathbf{c}^{-},\eta,\boldsymbol{\xi},\boldsymbol{\chi},\sigma\right)}$
(18)

where $g_{j}^{*}$, $s_{j}^{*}$ and $v_{j}^{*}$ represent the proposed values;
$g_{j}^{0}$, $s_{j}^{0}$ and $v_{j}^{0}$ are the old values; $\tilde{v}_{j+1}$
denotes $v_{j+1}$ under the proposed variable values (due to possible PYP
cluster label change or elimination). Specifically,
$\tilde{v}_{j+1}=(q^{*}_{g_{j+1}s_{j+1}}+1)$ if it belongs to no existing
cluster, where $q^{*}_{g_{j+1}s_{j+1}}$ is the number of clusters under the
proposed variable values and
$\boldsymbol{\phi}^{*}_{g_{j+1}s_{j+1}\tilde{v}_{j+1}}$ is the random effect
under the proposed variable values. Both the nominator and the denominator in
(18) can be calculated using (16) and (17), and by transposing $j$ to $(j+1)$.

##  9 Generalized form of the Sticky PYP

A Sticky PYP is defined by the following general properties:

  1. Property 1:

Set $\mathcal{G}$ comprises a countable number of generative groups. Each
group $g\in\mathcal{G}$ contains a countable number of group-specific regular
PYPs. The PYPs are identified by a combination of the group label $g$ and an
integer-valued state, $s\in\mathcal{S}$. In other words, the PYPs have
bivariate labels, $(g,s)\in\mathcal{G}\times\mathcal{S}$.

  2. Property 2:

For every $(g,s)\in\mathcal{G}\times\mathcal{S}$, let the corresponding PYP be
$\mathcal{W}_{gs}(d_{s},\alpha_{s},W_{s})$, of which the $T$-variate base
distribution $W_{s}$, discount parameter $d_{s}\in[0,1)$, and mass parameter
$\alpha_{s}>0$ depend on the state $s$. If set $\mathcal{S}$ contains multiple
states, assume that the base distributions $\\{W_{s}:s\in\mathcal{S}\\}$ are
such that two PYPs associated with unequal states almost surely have non-
intersecting sets of atoms.

Let distribution $\mathcal{P}_{gs}$ be an independent realizations of the
corresponding PYP:

$\mathcal{P}_{gs}\,\big{|}\,d_{s},\alpha_{s},W_{s}\stackrel{{\scriptstyle
indep}}{{\sim}}\mathcal{W}_{gs}(d_{s},\alpha_{s},W_{s}),\quad(g,s)\in\mathcal{G}\times\mathcal{S}.$

Consequently, for each $s$, the distributions
$\\{\mathcal{P}_{gs}:g\in\mathcal{G}\\}$ share the same atoms as base
distribution $W_{s}$; however, the probabilities associated with the atoms
depend on the group-state combination.

  3. Property 3:

For probe $j=1,\ldots,p$, the label of the PYP generating the random effect
$\boldsymbol{\theta}_{j}$ is denoted by
$(g_{j},s_{j})\in\mathcal{G}\times\mathcal{S}$, and
$\boldsymbol{\theta}_{j}\mid\mathcal{P}_{g_{j}s_{j}}\stackrel{{\scriptstyle
indep}}{{\sim}}\mathcal{P}_{g_{j}s_{j}}.$

  4. Property 4:

Given the group $g_{j}$ of the $j^{th}$ probe, the state has the distribution:

$s_{j}\mid g_{j}\sim\mathcal{Q}_{g_{j}}$

where, for every $g\in\mathcal{G}$, $\mathcal{Q}_{g}$ denotes a group-specific
probability mass function on the set $\mathcal{S}$; thus,
$\sum_{s\in\mathcal{S}}\mathcal{Q}_{g}(s)=1$.

  5. Property 5:

(Markov property) For the first probe, group $g_{1}$ follows a categorical
distribution, $\mathcal{F}_{1}$, on the set $\mathcal{G}$. For the subsequent
probes $j>1$, given the variables associated with the preceding probes, the
mass function of group variable $g_{j}$ generally depends on probe index $j$,
random vector $\boldsymbol{\theta}_{j-1}$, and inter-probe distance $e_{j-1}$:

$g_{j}\sim\mathcal{F}_{j\boldsymbol{\theta}_{j-1}e_{j-1}}$

where, for every $\boldsymbol{\theta}\in\mathcal{R}^{T}$ and $e>0$,
$\mathcal{F}_{j\boldsymbol{\theta}e}$ denotes a probability mass function on
the set $\mathcal{G}$, so that
$\sum_{g\in\mathcal{G}}\mathcal{F}_{j\boldsymbol{\theta}e}(g)=1$.

Some popular examples are presented in Table 2 of Supplementary Material. As
suggested by the form of distribution $\mathcal{F}_{j\boldsymbol{\theta}e}$ in
the table, rows 1–5 of Table 2 display zero order Sticky PYPs, e.g.,
hierarchical Pitman-Yor processes (HPYP) (Teh et al., 2006; Camerlenghi et
al., 2019). Rows 6 and 7 are first order Sticky PYPs involving multiple groups
and single states. The mass functions $\mathcal{F}_{j\boldsymbol{\theta}e}$ do
not depend on inter-probe distance $e$ in the examples.

For a detailed discussion of a specific mixture model, consider the
hierarchical Dirichlet process (HDP) of Teh et al. (2006), in which $p$
objects are organized into $K$ groups with known group memberships
$g_{1},\ldots,g_{p}$. Let global measure $W_{\infty}$ be a Dirichlet process
realization, and therefore, a countably infinite distribution in
$\mathcal{R}^{T}$. The random effects
$\boldsymbol{\theta}_{1},\ldots,\boldsymbol{\theta}_{p}$ are distributed as:

$\displaystyle\boldsymbol{\theta}_{j}\mid G_{g_{j}}$
$\displaystyle\stackrel{{\scriptstyle
indep}}{{\sim}}G_{g_{j}},\quad\text{$j=1,\ldots,p,$}$ $\displaystyle
G_{g}\mid\alpha,W_{\infty}$ $\displaystyle\stackrel{{\scriptstyle
i.i.d.}}{{\sim}}\text{DP}\bigl{(}\alpha,W_{\infty}\bigr{)},\quad\text{$g=1,\ldots,K$.}$
(19)

As summarized in the fourth row of Table 2, the model can be recast as a
generalized Sticky PYP as follows:

  1. Property 1:

The set of group labels, $\mathcal{G}=\\{1,\ldots,K\\}$. An HDP is a single-
state model with state set $\mathcal{S}=\\{1\\}$, i.e., $s_{j}=1$ for all $p$
objects. To simplify the notation, we henceforth drop the subscript $s$.

  2. Property 2:

For group $g=1,\ldots,K$, the associated PYP is
$\mathcal{W}_{g}\bigl{(}0,\alpha,W_{\infty}\bigr{)}=$
$\text{DP}\bigl{(}\alpha,W_{\infty}\bigr{)}$.

  3. Property 3:

Since distribution $\mathcal{P}_{g}$ in the generalized Sticky PYP formulation
is a random realization of
$\mathcal{W_{\infty}}_{g}\bigl{(}0,\alpha,W_{\infty}\bigr{)}=\text{DP}\bigl{(}\alpha,W_{\infty}\bigr{)}$,
we identify $\mathcal{P}_{g}=G_{g}$ in expression (19).

  4. Property 4:

An HDP is a single-state model; therefore, trivially,
$\mathcal{Q}_{g}=1_{\\{1\\}}$.

  5. Property 5:

Group variables $g_{1},\ldots,g_{p}$ are known. Hence,
$\mathcal{F}_{1}=1_{\\{g_{1}\\}}$ and
$\mathcal{F}_{j\boldsymbol{\theta},e}=1_{\\{g_{j}\\}}$, $j=2,\ldots,p$.

The 2R2CF used in differential methylation is a Sticky PYP with two groups and
two states. A key characteristic distinguishing it from the other first order
models, and indeed, all the models listed in Table 2, is that the 2R2CF makes
a probabilistic, rather than deterministic, assumption for group variable
$g_{j}$ that varies with the inter-probe distances. As noted in Section 2.1,
the 2R2CF behaves similarly to two-state hidden Markov models for very small
inter-probe distances and similarly to finite mixture models for relatively
large distances. For differential analysis, this offers a key advantage in
datasets with widely varying inter-probe distances by allowing the
differential and non-differential probes to have different cluster allocation
patterns depending on the differential statuses and distances of adjacent
probes. These interpretable features of the 2R2CF process are facilitated by
the two-group, two-state construct.

Model | $\mathcal{G}$ | $\mathcal{W}_{gs}(d_{s},\alpha_{s},W_{s})$ | $\mathcal{F}_{j\boldsymbol{\theta}e}$
---|---|---|---
Finite mixture | $\\{1\\}$ |  $\mathcal{W}(0,\alpha,W_{K})$, | $1_{\\{1\\}}$
Dirichlet process | $\\{1\\}$ | $\mathcal{W}(0,\alpha,W)$ | $1_{\\{1\\}}$
PYP | $\\{1\\}$ | $\mathcal{W}(d,\alpha,W)$ | $1_{\\{1\\}}$
HDP | $\mathcal{N}_{K}$ | $\mathcal{W}_{g}(0,\alpha,W_{\infty})$ |  $1_{\\{g_{j}\\}}$; known $g_{j}\in\mathcal{G}$
HPYP | $\mathcal{N}_{K}$ | $\mathcal{W}_{g}(d,\alpha,W_{\infty})$ |  $1_{\\{g_{j}\\}}$; known $g_{j}\in\mathcal{G}$
Finite HMM | $\mathcal{N}_{K}$ | $\mathcal{W}_{g}(0,\alpha,W_{K})$ |  Point mass at $\sum_{v=1}^{K}v\cdot\mathcal{I}(\boldsymbol{\theta}=\boldsymbol{\phi}_{v})$
HDP-HMM | $\mathcal{N}$ | $\mathcal{W}_{g}(0,\alpha,W_{\infty})$ |  Point mass at $\sum_{v=1}^{\infty}v\cdot\mathcal{I}(\boldsymbol{\theta}=\boldsymbol{\phi}_{v})$
Table 2: Examples of Sticky PYPs. The examples listed here correspond to a singleton state set, $\mathcal{S}=\\{1\\}$, and degenerate distribution, $\mathcal{Q}_{g}=1_{\\{1\\}}$. Set $\mathcal{N}$ represents the natural numbers, $\mathcal{N}_{K}=$ $\\{1,\ldots,K\\}$ for some $K>1$, $d\in[0,1)$, and $W$ is an arbitrary distribution in $\mathcal{R}^{T}$. Distribution $W_{K}$ is finite with $T$-variate atoms $\\{\boldsymbol{\phi}_{v}\\}_{v=1}^{K}$. Distribution $W_{\infty}$ is countably infinite with $T$-variate atoms $\\{\boldsymbol{\phi}_{v}\\}_{v=1}^{\infty}$. In an HDP or HDP-HMM, $W_{\infty}$ represents a Dirichlet process realization. For an HPYP, $W_{\infty}$ represents a PYP realization. Refer to Section 9 for the notation and further discussion.  |  | Low noise | High noise
---|---|---|---
|  | High correlation | No correlation | High correlation | No correlation
AUC | BayesDiff | 0.995 | 0.970 | 0.944 | 0.830
ANOVA | 0.964 | 0.949 | 0.888 | 0.864
Kruskal-Wallis | 0.958 | 0.944 | 0.878 | 0.854
COHCAP | 0.955 | 0.941 | 0.861 | 0.841
Methylkit | 0.955 | 0.942 | 0.869 | 0.848
BiSeq | 0.949 | 0.933 | 0.867 | 0.845
RADMeth | 0.959 | 0.944 | 0.876 | 0.852
$\textrm{AUC}_{20}$ | BayesDiff | 0.988 | 0.926 | 0.884 | 0.636
ANOVA | 0.891 | 0.851 | 0.684 | 0.628
Kruskal-Wallis | 0.876 | 0.833 | 0.657 | 0.597
COHCAP | 0.856 | 0.820 | 0.602 | 0.556
Methylkit | 0.858 | 0.819 | 0.633 | 0.579
BiSeq | 0.822 | 0.781 | 0.607 | 0.559
RADMeth | 0.871 | 0.833 | 0.639 | 0.588
$\textrm{AUC}_{10}$ | BayesDiff | 0.985 | 0.901 | 0.857 | 0.565
ANOVA | 0.849 | 0.804 | 0.586 | 0.529
Kruskal-Wallis | 0.827 | 0.773 | 0.551 | 0.493
COHCAP | 0.797 | 0.751 | 0.492 | 0.433
Methylkit | 0.794 | 0.744 | 0.523 | 0.458
BiSeq | 0.726 | 0.673 | 0.482 | 0.424
RADMeth | 0.817 | 0.775 | 0.531 | 0.471
Table 3: Areas under ROC curves for the different methods (rows) under the
four simulation scenarios (columns). See the text in the paper for further
discussion. Figure 6: Exploratory analysis of DNA methylation profiles of
upper GI cancer samples from TCGA. See the text in the paper for further
explanation.

(a) Low noise, high correlation

(b) Low noise, no correlation

(c) High noise, high correlation

(d) High noise, no correlation

Figure 7: ROC curves, averaged over 20 simulated datasets, for the seven
methods and the four simulation scenarios. The numbers in parentheses
represent AUCs.

(a) Histogram of the number of gene-specific CpG sites

(b) 95% credible intervals for lower bounds of log-Bayes factors of first
order versus zero-order models. The intervals whose lower limits are positive
are marked in red.

Figure 8: Data analysis plots. See the Section 5 text for further explanation.
BayesDiff
---
Gene | DM | Largest | Gene | DM | Largest
Proportion | Difference | Proportion | Difference
EDNRB | 1.00 | STAD$\uparrow$ LIHC$\downarrow$ | SLITRK1 | 0.96 | STAD$\uparrow$ LIHC$\downarrow$
PCLO | 1.00 | PAAD$\uparrow$ LIHC$\downarrow$ | FLRT2 | 0.96 | STAD$\uparrow$ LIHC$\downarrow$
PREX2 | 1.00 | STAD$\uparrow$ LIHC$\downarrow$ | KCNA1 | 0.96 | STAD$\uparrow$ LIHC$\downarrow$
SLIT2 | 1.00 | STAD$\uparrow$ LIHC$\downarrow$ | TRPA1 | 0.96 | STAD$\uparrow$ LIHC$\downarrow$
SLITRK2 | 1.00 | STAD$\uparrow$ LIHC$\downarrow$ | ADCY8 | 0.96 | STAD$\uparrow$ LIHC$\downarrow$
SORCS3 | 1.00 | STAD$\uparrow$ LIHC$\downarrow$ | CTNNA2 | 0.95 | PAAD$\uparrow$ LIHC$\downarrow$
SPHKAP | 1.00 | ESCA$\uparrow$ LIHC$\downarrow$ | GRIA2 | 0.95 | STAD$\uparrow$ LIHC$\downarrow$
SPTA1 | 1.00 | PAAD$\uparrow$ LIHC$\downarrow$ | ADGRL3 | 0.94 | PAAD$\uparrow$ LIHC$\downarrow$
UNC13C | 1.00 | PAAD$\uparrow$ LIHC$\downarrow$ | LRRC7 | 0.94 | STAD$\uparrow$ LIHC$\downarrow$
XIRP2 | 1.00 | PAAD$\uparrow$ LIHC$\downarrow$ | ERBB4 | 0.93 | PAAD$\uparrow$ LIHC$\downarrow$
ZNF804B | 1.00 | STAD$\uparrow$ LIHC$\downarrow$ | PCDH10 | 0.93 | STAD$\uparrow$ LIHC$\downarrow$
TSHZ3 | 0.97 | STAD$\uparrow$ LIHC$\downarrow$ | SOX11 | 0.93 | STAD$\uparrow$ LIHC$\downarrow$
MYO3A | 0.97 | STAD$\uparrow$ LIHC$\downarrow$ | NLGN4X | 0.93 | PAAD$\uparrow$ LIHC$\downarrow$
ABCC9 | 0.97 | ESCA$\uparrow$ LIHC$\downarrow$ | NBEA | 0.93 | PAAD$\uparrow$ LIHC$\downarrow$
EPB41L3 | 0.97 | STAD$\uparrow$ LIHC$\downarrow$ | CNTN1 | 0.92 | STAD$\uparrow$ LIHC$\downarrow$
FBN2 | 0.97 | STAD$\uparrow$ LIHC$\downarrow$ | GRM5 | 0.92 | PAAD$\uparrow$ LIHC$\downarrow$
PCDH17 | 0.96 | STAD$\uparrow$ LIHC$\downarrow$ | PTPRZ1 | 0.91 | STAD$\uparrow$ PAAD$\downarrow$
CDH8 | 0.96 | STAD$\uparrow$ LIHC$\downarrow$ | EPHA5 | 0.91 | STAD$\uparrow$ LIHC$\downarrow$
Table 4: For BayesDiff, genes with the overall proportion of differentially
methylated CpG sites exceeding 0.9. The “Largest Difference” column displays
which pair-wise difference between the four cancer types is the largest, with
the symbol “$\uparrow$ ($\downarrow$)” indicating higher (lower) methylation
level for one cancer type relative to the other. BayesDiff
---
Gene | Number of probes | Gene | Number of probes
ABCC9 | 36 | NLGN4X | 28
ADCY8 | 24 | PCDH10 | 30
ADGRL3 | 35 | PCDH17 | 27
CDH8 | 26 | PCLO | 20
CNTN1 | 25 | PREX2 | 21
CTNNA2 | 85 | PTPRZ1 | 23
EDNRB | 56 | SLIT2 | 29
EPB41L3 | 34 | SLITRK1 | 26
EPHA5 | 22 | SLITRK2 | 28
ERBB4 | 30 | SORCS3 | 39
FBN2 | 30 | SOX11 | 44
FLRT2 | 51 | SPHKAP | 17
GRIA2 | 20 | SPTA1 | 16
GRM5 | 25 | TRPA1 | 25
KCNA1 | 25 | TSHZ3 | 39
LRRC7 | 34 | UNC13C | 11
MYO3A | 38 | XIRP2 | 19
NBEA | 55 | ZNF804B | 21
Table 5: Number of included CpG sites for the top methylated genes listed in
Table 4 for BayesDiff. BayesDiff and ANOVA
---
Gene | DM Proportion | Largest Difference | DM Proportion | Largest Difference
BayesDiff | BayesDiff | ANOVA | ANOVA
PCLO | 1.00 | PAAD↑ LIHC↓ | 1.00 | PAAD↑ LIHC↓
UNC13C | 1.00 | PAAD↑ LIHC↓ | 1.00 | PAAD↑ ESCA↓
XIRP2 | 1.00 | PAAD↑ LIHC↓ | 1.00 | PAAD↑ STAD↓
ZNF804B | 1.00 | STAD↑ LIHC↓ | 1.00 | ESCA↑ PAAD↓
FBN2 | 0.97 | STAD↑ LIHC↓ | 0.93 | PAAD↑ STAD↓
FLRT2 | 0.96 | STAD↑ LIHC↓ | 0.92 | STAD↑ PAAD↓
ERBB4 | 0.93 | PAAD↑ LIHC↓ | 0.93 | STAD↑ PAAD↓
PTPRZ1 | 0.91 | STAD↑ PAAD↓ | 0.91 | PAAD↑ LIHC↓
Table 6: For both BayesDiff and ANOVA methods, the common set of genes
detected with the overall proportion of differentially methylated CpG sites
exceeding 0.9. For each method, the “Largest Difference” column displays which
pair-wise difference between the four cancer types is the largest, with the
symbol “$\uparrow$ ($\downarrow$)” indicating higher (lower) methylation level
for one cancer type relative to the other. ANOVA, not BayesDiff
---
Gene | DM | Largest | Gene | DM | Largest
Proportion | Difference | Proportion | Difference
LRP12 | 1.00 | LIHC↑ STAD↓ | EPHA6 | 0.95 | STAD↑ PAAD↓
PTPRM | 1.00 | STAD↑ LIHC↓ | RELN | 0.95 | STAD↑ LIHC↓
SALL1 | 0.96 | STAD↑ LIHC↓ | COL6A6 | 0.94 | STAD↑ LIHC↓
TIAM1 | 0.96 | LIHC↑ PAAD↓ | MUC16 | 0.94 | PAAD↑ LIHC↓
MXRA5 | 0.96 | LIHC↑ PAAD↓ | HCN1 | 0.91 | PAAD↑ ESCA↓
ZFHX4 | 0.96 | LIHC↑ STAD↓ | RIMS2 | 0.91 | STAD↑ LIHC↓
Table 7: For only ANOVA but not BayesDiff, genes detected with the overall
proportion of differentially methylated CpG sites exceeding 0.9. The “Largest
Difference” column displays which pair-wise difference between the four cancer
types is the largest, with the symbol “$\uparrow$ ($\downarrow$)” indicating
higher (lower) methylation level for one cancer type relative to the other.
BayesDiff, not ANOVA
---
Gene | DM | Largest | Gene | DM | Largest
Proportion | Difference | Proportion | Difference
EDNRB | 1.00 | STAD↑ LIHC↓ | KCNA1 | 0.96 | STAD↑ LIHC↓
PREX2 | 1.00 | STAD↑ LIHC↓ | TRPA1 | 0.96 | STAD↑ LIHC↓
SLIT2 | 1.00 | STAD↑ LIHC↓ | ADCY8 | 0.96 | STAD↑ LIHC↓
SLITRK2 | 1.00 | STAD↑ LIHC↓ | CTNNA2 | 0.95 | PAAD↑ LIHC↓
SORCS3 | 1.00 | STAD↑ LIHC↓ | GRIA2 | 0.95 | STAD↑ LIHC↓
SPHKAP | 1.00 | ESCA↑ LIHC↓ | ADGRL3 | 0.94 | PAAD↑ LIHC↓
SPTA1 | 1.00 | PAAD↑ LIHC↓ | LRRC7 | 0.94 | STAD↑ LIHC↓
TSHZ3 | 0.97 | STAD↑ LIHC↓ | PCDH10 | 0.93 | STAD↑ LIHC↓
MYO3A | 0.97 | STAD↑ LIHC↓ | SOX11 | 0.93 | STAD↑ LIHC↓
ABCC9 | 0.97 | ESCA↑ LIHC↓ | NLGN4X | 0.93 | PAAD↑ LIHC↓
EPB41L3 | 0.97 | STAD↑ LIHC↓ | NBEA | 0.93 | PAAD↑ LIHC↓
PCDH17 | 0.96 | STAD↑ LIHC↓ | CNTN1 | 0.92 | STAD↑ LIHC↓
CDH8 | 0.96 | STAD↑ LIHC↓ | GRM5 | 0.92 | PAAD↑ LIHC↓
SLITRK1 | 0.96 | STAD↑ LIHC↓ | EPHA5 | 0.91 | STAD↑ LIHC↓
Table 8: For only BayesDiff but not ANOVA, genes detected with overall
proportion of differentially methylated CpG sites exceeding 0.9. The “Largest
Difference” column displays which pair-wise difference between the four cancer
types is the largest, with the symbol “$\uparrow$ ($\downarrow$)” indicating
higher (lower) methylation level for one cancer type relative to the other.

(a) Gene TP53

(b) Gene TTN

Figure 9: Detailed differential methylation results for the top 2 mutated
genes. For each gene, the upper panel shows the mean methylation levels. The
middle panel shows the posterior probabilities of each CpG site being
differentially methylated, with solid points representing differential
methylation and dashed line denoting the corresponding cutoff value. The lower
panel indicates the largest pairwise difference between the 4 cancer types.
Symbols 1–4 in the lower panel represent GI cancer types STAD, LIHC, ESCA and
PAAD, respectively. The vertical dotted lines represent the gene boundaries.
The arrow at the top indicates the transcription direction.

(a) First moments for gene TP53

(b) Second moments for gene TP53

(c) First moments for gene TTN

(d) Second moments for gene TTN

Figure 10: Comparison of predicted and actual sample moments for the 2 top
mutated genes. ptRestaurant and cuisine of customer $j$Cuisine menusptSubject-
specific random
effects$\eta$$\alpha_{s}$$d_{s}$$\beta$$e_{j-1}$$r_{j}$$v_{j}$$W_{s}$$G$$\rho_{2}$$g_{j}$$\boldsymbol{\phi}_{gsv}$$\mu_{G},\tau_{G}^{2}$$\gamma$$s_{j}$$\boldsymbol{\theta}_{j}$$\mathbf{x}_{j}$$\mathbf{z}_{j}$$\xi_{i}$$\epsilon_{i}$$\chi_{j}$$\sigma^{2}$$b_{i}$$\tau_{\epsilon}^{2}$
Figure 11: Directed acyclic graph of the BayesDiff model showing the
relationships between the model parameters. Solid rectangles represent the
data and model parameters. Dashed rectangles represent predetermined
constants. Solid arrows represent stochastic relationships, and dashed arrows
represent deterministic relationships.

## References

  * (1)
  * Akalin et al. (2012) Akalin, A., Kormaksson, M., Li, S., Garrett-Bakelman, F. E., Figueroa, M. E., Melnick, A. & Mason, C. E. (2012), ‘methylkit: a comprehensive r package for the analysis of genome-wide dna methylation profiles’, Genome Biology 13(10), R87.
  * Bae et al. (2017) Bae, H., Kim, B., Lee, H., Lee, S., Kang, H.-S. & Kim, S. J. (2017), ‘Epigenetically regulated fibronectin leucine rich transmembrane protein 2 (flrt2) shows tumor suppressor activity in breast cancer cells’, Scientific Reports 7(1), 272.
  * Baker et al. (2008) Baker, T. A., Bell, S. P., Gann, A., Levine, M., Losick, R. & Inglis, C. (2008), Molecular biology of the gene, San Francisco, CA, USA:: Pearson/Benjamin Cummings.
  * Basu & Chib (2003) Basu, S. & Chib, S. (2003), ‘Marginal likelihood and bayes factors for dirichlet process mixture models’, Journal of the American Statistical Association 98(461), 224–235.
  * Benjamini & Hochberg (1995) Benjamini, Y. & Hochberg, Y. (1995), ‘Controlling the false discovery rate: a practical and powerful approach to multiple testing’, Journal of the Royal Statistical Society: Series B (Methodological) pp. 289–300.
  * Blei & Jordan (2005) Blei, D. M. & Jordan, M. I. (2005), ‘Variational inference for dirichlet process mixtures’, Bayesian Analysis 1, 1–23.
  * Camerlenghi et al. (2019) Camerlenghi, F., Lijoi, A., Orbanz, P. & Prünster, I. (2019), ‘Distribution theory for hierarchical processes’, The Annals of Statistics 47(1), 67 – 92\.
  * Chib (1995) Chib, S. (1995), ‘Marginal likelihood from the gibbs output’, Journal of the American Statistical Association 90(432), 1313–1321.
  * Chib & Winkelmann (2001) Chib, S. & Winkelmann, R. (2001), ‘Markov chain monte carlo analysis of correlated count data’, Journal of Business & Economic Statistics 19(4), 428–435.
  * De Blasi et al. (2013) De Blasi, P., Favaro, S., Lijoi, A., Mena, R. H., Prünster, I. & Ruggiero, M. (2013), ‘Are gibbs-type priors the most natural generalization of the dirichlet process?’, IEEE transactions on pattern analysis and machine intelligence 37(2), 212–229.
  * Dolzhenko & Smith (2014) Dolzhenko, E. & Smith, A. D. (2014), ‘Using beta-binomial regression for high-precision differential methylation analysis in multifactor whole-genome bisulfite sequencing experiments’, BMC bioinformatics 15(1), 215.
  * Du et al. (2010) Du, P., Zhang, X., Huang, C.-C., Jafari, N., Kibbe, W. A., Hou, L. & Lin, S. M. (2010), ‘Comparison of beta-value and m-value methods for quantifying methylation levels by microarray analysis’, BMC bioinformatics 11, 1–9.
  * Dunson et al. (2008) Dunson, D. B., Herring, A. H. & Engel, S. M. (2008), ‘Bayesian selection and clustering of polymorphisms in functionally-related genes’, Journal of the American Statistical Association 103, 534–546.
  * Dunson & Park (2008) Dunson, D. B. & Park, J.-H. (2008), ‘Kernel stick-breaking processes’, Biometrika 95, 307–323.
  * Durbin et al. (2002) Durbin, B., Hardin, J., Hawkins, D. & Rocke, D. (2002), ‘A variance-stabilizing transformation for gene-expression microarray data’, Bioinformatics 18, S105––S110.
  * Eckhardt et al. (2006) Eckhardt, F., Lewin, J., Cortese, R., Rakyan, V. K., Attwood, J., Burger, M., Burton, J., Cox, T. V., Davies, R., Down, T. A. et al. (2006), ‘Dna methylation profiling of human chromosomes 6, 20 and 22’, Nature Genetics 38(12), 1378.
  * Escobar (1994) Escobar, M. D. (1994), ‘Estimating normal means with a dirichlet process prior’, Journal of the American Statistical Association 89(425), 268–277.
  * Escobar & West (1995) Escobar, M. D. & West, M. (1995), ‘Bayesian density estimation and inference using mixtures’, Journal of the American statistical association 90(430), 577–588.
  * Feinberg & Tycko (2004) Feinberg, A. P. & Tycko, B. (2004), ‘The history of cancer epigenetics’, Nature reviews. Cancer 4(2), 143.
  * Feng et al. (2014) Feng, H., Conneely, K. N. & Wu, H. (2014), ‘A bayesian hierarchical model to detect differentially methylated loci from single nucleotide resolution sequencing data’, Nucleic Acids Research 42(8), e69–e69.
  * Ferguson (1973a) Ferguson, T. S. (1973a), ‘A bayesian analysis of some nonparametric problems’, The Annals of Statistics pp. 209–230.
  * Ferguson (1973b) Ferguson, T. S. (1973b), ‘A bayesian analysis of some nonparametric problems’, Annals of Statistics 1, 209–223.
  * Fox et al. (2011) Fox, E., Sudderth, E., Jordan, M. & Willsky, A. (2011), ‘The sticky hdp-hmm: Bayesian nonparametric hidden markov models with persistent states’, Annals of Applied Statistics 5, 1020–1056.
  * Frühwirth-Schnatter (2006) Frühwirth-Schnatter, S. (2006), Finite Mixture and Markov Switching Models, New York: Springer.
  * Gnedin & Pitman (2006) Gnedin, A. & Pitman, J. (2006), ‘Exchangeable gibbs partitions and stirling triangles’, Journal of Mathematical Sciences 138, 5674–5685.
  * Grossman et al. (2016) Grossman, R. L., Heath, A. P., Ferretti, V., Varmus, H. E., Lowy, D. R., Kibbe, W. A. & Staudt, L. M. (2016), ‘Toward a shared vision for cancer genomic data’, New England Journal of Medicine 375(12), 1109–1112.
  * Guha (2010) Guha, S. (2010), ‘Posterior simulation in countable mixture models for large datasets’, Journal of the American Statistical Association 105(490), 775–786.
  * Guha & Baladandayuthapani (2016) Guha, S. & Baladandayuthapani, V. (2016), ‘A nonparametric bayesian technique for high-dimensional regression’, Electronic Journal of Statistics 10, 3374–3424.
  * Hansen et al. (2012) Hansen, K. D., Langmead, B. & Irizarry, R. A. (2012), ‘Bsmooth: from whole genome bisulfite sequencing reads to differentially methylated regions’, Genome Biology 13(10), R83.
  * Hebestreit et al. (2013) Hebestreit, K., Dugas, M. & Klein, H.-U. (2013), ‘Detection of significantly differentially methylated regions in targeted bisulfite sequencing data’, Bioinformatics 29(13), 1647–1653.
  * Irizarry et al. (2008) Irizarry, R. A., Ladd-Acosta, C., Carvalho, B., Wu, H., Brandenburg, S. A., Jeddeloh, J. A., Wen, B. & Feinberg, A. P. (2008), ‘Comprehensive high-throughput arrays for relative methylation (charm)’, Genome Research 18(5), 780–790.
  * Irizarry et al. (2009) Irizarry, R. A., Ladd-Acosta, C., Wen, B., Wu, Z., Montano, C., Onyango, P., Cui, H., Gabo, K., Rongione, M., Webster, M. et al. (2009), ‘Genome-wide methylation analysis of human colon cancer reveals similar hypo-and hypermethylation at conserved tissue-specific cpg island shores’, Nature Genetics 41(2), 178.
  * Ishwaran & James (2003) Ishwaran, H. & James, L. F. (2003), ‘Generalized weighted chinese restaurant processes for species sampling mixture models’, Statistica Sinica 13, 1211–1235.
  * Ishwaran & Zarepour (2002) Ishwaran, H. & Zarepour, M. (2002), ‘Dirichlet prior sieves in finite normal mixtures’, Statistica Sinica 12, 941–963.
  * Jaffe et al. (2012) Jaffe, A. E., Murakami, P., Lee, H., Leek, J. T., Fallin, M. D., Feinberg, A. P. & Irizarry, R. A. (2012), ‘Bump hunting to identify differentially methylated regions in epigenetic epidemiology studies’, International Journal of Epidemiology 41(1), 200–209.
  * Kim et al. (2006) Kim, S., Tadesse, M. G. & Vannucci, M. (2006), ‘Variable selection in clustering via dirichlet process mixture models’, Biometrika 93, 877–893.
  * Leek et al. (2010) Leek, J. T., Scharpf, R. B., Bravo, H. C., Simcha, D., Langmead, B., Johnson, W. E., Geman, D., Baggerly, K. & Irizarry, R. A. (2010), ‘Tackling the widespread and critical impact of batch effects in high-throughput data’, Nature Reviews Genetics 11(10).
  * Li et al. (2016) Li, D., Wang, X., Lin, L. & Dey, D. K. (2016), ‘Flexible link functions in nonparametric binary regression with gaussian process priors’, Biometrics 72(3), 707––719.
  * Lijoi et al. (2007a) Lijoi, A., Mena, R. & Prünster, I. (2007a), ‘Bayesian nonparametric estimation of the probability of discovering new species’, Biometrika 94, 769–786.
  * Lijoi et al. (2007b) Lijoi, A., Mena, R. & Prünster, I. (2007b), ‘Controlling the reinforcement in bayesian nonparametric mixture models’, Journal of the Royal Statistical Society: Series B (Statistical Methodology) 69, 715–740.
  * Lijoi & Prünster (2010) Lijoi, A. & Prünster, I. (2010), Models beyond the Dirichlet process, Cambridge Series in Statistical and Probabilistic Mathematics, pp. 80–136.
  * MacEachern (1994) MacEachern, S. N. (1994), ‘Estimating normal means with a conjugate style dirichlet process prior’, Communications in Statistics-Simulation and Computation 23(3), 727–741.
  * Medvedovic et al. (2004) Medvedovic, M., Yeung, K. Y. & Bumgarner, R. E. (2004), ‘Bayesian mixture model based clustering of replicated microarray data’, Bioinformatics 20, 1222–1232.
  * Müller & Mitra (2013) Müller, P. & Mitra, R. (2013), ‘Bayesian nonparametric inference–why and how’, Bayesian Analysis (Online) 8(2).
  * Murphy et al. (2013) Murphy, S. J., Hart, S. N., Lima, J. F., Kipp, B. R., Klebig, M., Winters, J. L., Szabo, C., Zhang, L., Eckloff, B. W., Petersen, G. M. et al. (2013), ‘Genetic alterations associated with progression from pancreatic intraepithelial neoplasia to invasive pancreatic tumor’, Gastroenterology 145(5), 1098–1109.
  * Newton et al. (2004) Newton, M. A., Noueiry, A., Sarkar, D. & Ahlquist, P. (2004), ‘Detecting differential gene expression with a semiparametric hierarchical mixture method’, Biostatistics 5(2), 155–176.
  * Park et al. (2014) Park, Y., Figueroa, M. E., Rozek, L. S. & Sartor, M. A. (2014), ‘Methylsig: a whole genome dna methylation analysis pipeline’, Bioinformatics 30(17), 2414–2422.
  * Perman et al. (1992) Perman, M., Pitman, J. & Yor, M. (1992), ‘Size-biased sampling of poisson point processes and excursions’, Probability Theory and Related Fields 92(1), 21–39.
  * Rackham et al. (2015) Rackham, O. J., Dellaportas, P., Petretto, E. & Bottolo, L. (2015), ‘Wgbssuite: simulating whole-genome bisulphite sequencing data and benchmarking differential dna methylation analysis tools’, Bioinformatics 31(14), 2371–2373.
  * Rodriguez et al. (2008) Rodriguez, A., B., D. D. & Gelfand, A. E. (2008), ‘The nested dirichlet process (with discussion)’, Journal of the American Statistical Association 103, 1131–1144.
  * Saito et al. (2014) Saito, Y., Tsuji, J. & Mituyama, T. (2014), ‘Bisulfighter: accurate detection of methylated cytosines and differentially methylated regions’, Nucleic Acids Research p. gkt1373.
  * Sethuraman (1994) Sethuraman, J. (1994), ‘A constructive definition of dirichlet priors’, Statistica Sinica pp. 639–650.
  * Siegel et al. (2017) Siegel, R. L., Miller, K. D. & Jemal, A. (2017), ‘Cancer statistics, 2017’, CA: A Cancer Journal for Clinicians 67(1), 7–30.
  * Song et al. (2013) Song, Q., Decato, B., Hong, E. E., Zhou, M., Fang, F., Qu, J., Garvin, T., Kessler, M., Zhou, J. & Smith, A. D. (2013), ‘A reference methylome database and analysis pipeline to facilitate integrative and comparative epigenomics’, PloS One 8(12), e81148.
  * Subramaniam & Hsiao (2012) Subramaniam, S. & Hsiao, G. (2012), ‘Gene-expression measurement: variance-modeling considerations for robust data analysis’, Nature Immunology 13(3), 199–203.
  * Sun et al. (2014) Sun, D., Xi, Y., Rodriguez, B., Park, H. J., Tong, P., Meong, M., Goodell, M. A. & Li, W. (2014), ‘Moabs: model based analysis of bisulfite sequencing data’, Genome Biology 15(2), R38.
  * Tao et al. (2012) Tao, K., Wu, C., Wu, K., Li, W., Han, G., Shuai, X. & Wang, G. (2012), ‘Quantitative analysis of promoter methylation of the ednrb gene in gastric cancer’, Medical Oncology 29(1), 107–112.
  * Teh et al. (2006) Teh, Y. W., Jordan, M. I., Beal, M. J. & Blei, D. M. (2006), ‘Hierarchical dirichlet processes’, Journal of the American Statistical Association 101, 1566–1581.
  * Tomlinson & Escobar (2003) Tomlinson, G. & Escobar, M. (2003), ‘Analysis of densities’, Talk given at the Joint Statistical Meeting 103, 1131–1144.
  * Tsunoda et al. (2009) Tsunoda, S., Smith, E., De Young, N. J., Wang, X., Tian, Z.-Q., Liu, J.-F., Jamieson, G. G. & Drew, P. A. (2009), ‘Methylation of cldn6, fbn2, rbp1, rbp4, tfpi2, and tmeff2 in esophageal squamous cell carcinoma’, Oncology Reports 21(4), 1067–1073.
  * Vedeld et al. (2017) Vedeld, H. M., Goel, A. & Lind, G. E. (2017), Epigenetic biomarkers in gastrointestinal cancers: The current state and clinical perspectives, in ‘Seminars in cancer biology’, Elsevier.
  * Vrana et al. (2018) Vrana, D., Hlavac, V., Brynychova, V., Vaclavikova, R., Neoral, C., Vrba, J., Aujesky, R., Matzenauer, M., Melichar, B. & Soucek, P. (2018), ‘Abc transporters and their role in the neoadjuvant treatment of esophageal cancer’, International Journal of Molecular Sciences 19(3), 868.
  * Wang et al. (2012) Wang, D., Yan, L., Hu, Q., Sucheston, L. E., Higgins, M. J., Ambrosone, C. B., Johnson, C. S., Smiraglia, D. J. & Liu, S. (2012), ‘Ima: an r package for high-throughput analysis of illumina’s 450k infinium methylation data’, Bioinformatics 28(5), 729–730.
  * Warden et al. (2013) Warden, C. D., Lee, H., Tompkins, J. D., Li, X., Wang, C., Riggs, A. D., Yu, H., Jove, R. & Yuan, Y.-C. (2013), ‘Cohcap: an integrative genomic pipeline for single-nucleotide resolution dna methylation analysis’, Nucleic Acids Research 41(11), e117–e117.
  * Yu & Sun (2016) Yu, X. & Sun, S. (2016), ‘Hmm-dm: identifying differentially methylated regions using a hidden markov model’, Statistical Applications in Genetics and Molecular Biology 15(1), 69–81.
  * Zeger & Karim (1991) Zeger, S. L. & Karim, M. R. (1991), ‘Generalized linear models with random effects: A gibbs sampling approach’, Journal of the American Statistical Association 86, 79–86.