11institutetext: Columbia University
11email<EMAIL_ADDRESS>
# Evolving Interpretable Visual Classifiers with Large Language Models

Mia Chiquier  Utkarsh Mall  Carl Vondrick

###### Abstract

Multimodal pre-trained models, such as CLIP, are popular for zero-shot
classification due to their open-vocabulary flexibility and high performance.
However, vision-language models, which compute similarity scores between
images and class labels, are largely black-box, with limited interpretability,
risk for bias, and inability to discover new visual concepts not written down.
Moreover, in practical settings, the vocabulary for class names and attributes
of specialized concepts will not be known, preventing these methods from
performing well on images uncommon in large-scale vision-language datasets. To
address these limitations, we present a novel method that discovers
interpretable yet discriminative sets of attributes for visual recognition. We
introduce an evolutionary search algorithm that uses a large language model
and its in-context learning abilities to iteratively mutate a concept
bottleneck of attributes for classification. Our method produces state-of-the-
art, interpretable fine-grained classifiers. We outperform the latest
baselines by $18.4\%$ on five fine-grained iNaturalist datasets and by
$22.2\%$ on two KikiBouba datasets, despite the baselines having access to
privileged information about class names.

###### Keywords:

Visual Recognition, Interpretable Representations

##  1 Introduction

Multimodal foundation models like CLIP [1] obtain excellent performance on
many visual recognition tasks due to their flexibility to represent open-
vocabulary classes. These models have the potential to impact many scientific
applications, where computer vision systems could automate recognition in
specialized domains. However, since foundation models are neural networks,
they are largely black-box and we therefore have no means to explain or audit
the predictions they produce, limiting their trust. Moreover, given that
foundation models are trained on large corpora of web-scraped data [1, 2],
they are not optimized to represent rare and fine-grained concepts, such as
images from various scientific fields, which are infrequently described on the
internet.

The computer vision community has been building interpretable models by
integrating language, where classifiers are constructed with a bottleneck of
sparse or discrete attributes [3, 4, 5, 6, 7]. Language-based approaches have
the benefit of being interpretable. However, they rely on attributes that are
either hand-designed, requiring expert knowledge, or extracted from external
sources, such as large language models (LLMs). Since the attributes are not
learned, they often obtain poor performance on specialized classes
infrequently discussed in web-scale training sets.

Figure 1: Learning Interpretable Classifiers. Can we find text attributes for
a concept by looking at the images without their class names? LLM-Mutate is a
framework that learns sets of maximally discriminative visual attributes per
class without access to class names or any form of prior knowledge.

In this paper, we propose a framework to learn interpretable visual
representations from images. Our approach can discover discrete,
discriminative attributes from purely visual training sets of specialized
concepts, and even unwritten concepts that do not appear on the internet.
Solving this problem has typically been challenging because the loss is not
differentiable with respect to discrete attributes, limiting the application
of modern deep learning methods. Gradient-free methods, such as evolutionary
search [8], are able to optimize over discrete search spaces, but they have
not scaled to large problems in the past.

Our approach integrates large language models with evolutionary algorithms in
order to learn discriminative and interpretable attributes for visual
recognition. Evolutionary algorithms work by maintaining a set of candidates,
randomly mutating the candidates, and discarding the poorly performing
solutions according to an objective function. The mutation step is the
bottleneck for large-scale problems, such as in object recognition, because
the search space is large and the lack of gradient information means the
mutations are not guaranteed to drive the optimization towards rapid
convergence.

We overcome these bottlenecks by replacing the mutation step with a large
language model instead, whose in-context learning abilities are able to find
patterns in-between candidates and predict strong mutations that reduce the
objective. We evaluate our method on images from specialized scientific domain
that have been infrequently discussed on the internet due to the nature of
their specificity, with the iNaturalist dataset [9]. We chose five families of
plants and animals within iNaturalist, each containing between five and six
species, and evaluated the fine-grained classification performance on each
family. We outperform all baselines, on average by $18.4\%$ per family
dataset. We also evaluated our method of imaginary concepts that do not exist,
and consequently have been hardly discussed in training sets before. Following
the KikiBouba experiments [10] where people associate imaginary objects with
non-existent words, we learn interpretable attributes that achieve strong
discriminative performance, outperforming baselines by an average of $22.2\%$.

Our primary contribution is a framework for learning interpretable visual
recognition systems from image data. We propose to tightly integrate
evolutionary search with large language models, allowing us to efficiently
learn discrete, discriminative attributes that are interpretable. The
remainder of the paper will discuss the related work in this area, describe
the method, and present experimental results on multiple datasets. Our code,
data, and models are available at https://llm-mutate.cs.columbia.edu.

##  2 Related Work

We briefly review related work in interpretable methods, vision-language
models and evolutionary search.

Concept Bottlenecks. A common approach for interpretability is to create a
bottleneck in the classifier, where interpretable attributes are first
predicted before classifying the object category [3, 4, 6, 11, 12, 13].
Concept bottleneck models have been comprehensively studied in the domain of
zero-shot learning [14, 15, 16, 17, 18]. However, all these methods require a
significant amount of annotations for the intermediate attributes to perform
well while being interpretable. Our method instead finds interpretable
intermediate concepts by learning from visual data.

Post-hoc Interpretabililty. Several post-hoc methods have been proposed to
interpret deep models as well. Methods such as GradCAM [19] rely on the
activation maps to provide explanations for the classification on an image
[20, 21, 22, 23, 24, 25]. Generating counterfactual examples is another way to
interpret a model [26, 27, 28, 29]. Other work has looked at sample importance
for explanations [30, 31, 32, 33, 34, 35, 36]. All these methods look at
explanations from the perspective of a single image. In contrast to these
approaches, several methods aim at understanding the behaviors of individual
neurons in the trained model [37, 38, 39, 40, 41]. One major limitation of
post-hoc interpretability methods they can only explain the receptive field of
the model, which requires some interpretation on the user’s part. Concept or
attribute bottleneck models on the other hand can justify the classification
by providing scores for individual attributes.

Vision-Language Models. Large-scale vision-language models (VLMs), such as
CLIP [1], have bridged the visual and language modalities through contrastive
learning. VLMs estimate the similarity between text and an image, leading to
many downstream applications, such as zero-shot classification with open-ended
language. Several improvements have been proposed to the VLMs, such as
training with noisy data [42, 43, 44], better training strategies [45, 46, 47,
48, 49, 50, 51, 52], concept localization [53], multiple modalities [54],
grounding abilities [55, 56], or generative capabilities [57, 58, 59, 60].
Another advantage of VLMs is that classification or retrieval can be done by
using a list of attribute descriptors in conjunction with the class names [5,
61, 62]. Classification by using descriptions results in an intermediate step
of interpretability, as these lists of descriptions are essentially concept
bottlenecks. Waffle-CLIP [7] proposed using random words along with these
concepts leads to more robust concept embedding. However, all these methods
have access to a knowledge base (large language models in this case), which
they use to retrieve descriptions for categories using the category name. Such
approaches fail if the categories are never discussed or are infrequently
mentioned on the internet. Since our approach does not rely on the name of the
category, we perform better on specialized and esoteric categories.

LLMs as Optimizers. Using a large language model to generate mutations with
evolutionary search has been implemented for the task of code generation [63].
We follow a similar approach, except for the task of visual recognition. In
addition, their fitness functions are a series of non-learning-based metrics.
We instead use another foundation model, a vision-language model, as the
fitness function, thereby scaling to open-ended visual concepts. Another
instance in which LLM’s are used as optimizers is for prompt discovery [64].
Similarly, concurrent work leverages LLMs as optimizers to find attributes to
improve visual classification, but crucially, they provide the class name as
an input to the large language model [65, 66]. As such, they use privileged
information that we do not assume.

##  3 Discovering Visual Classifiers

Figure 2: Method. LLM-Mutate is an evolutionary algorithm that learns sets of
discrete language attributes per class. The mutation and cross-over
operations, which are mechanisms to introduce new parameter hypotheses, are
replaced by a large language model that uses in-context learning over past
attributes and their scores to iteratively generate better attributes.

###  3.1 Model

Given an image $x$, our goal is to predict its category label $y$. We create a
concept bottleneck model $f_{c}(x)$ that scores whether the image $x$ contains
category $c$. Once optimized, the concept bottleneck model should produce a
high score for $c=y$ and a low score otherwise. To score the class $c$ from
the image $x$, we average over a bottleneck of discrete attributes:

$\displaystyle
f_{c}(x;\mathcal{D})=\frac{1}{\left|\mathcal{D}(c)\right|}\sum_{d_{i}\in\mathcal{D}(c)}\phi(d_{i},x)$
(1)

where $\mathcal{D}(c)$ is the set of attributes to recognize class $c$, and
$\phi(d_{i},x)$ is the score from a vision-language model (such as CLIP) to
detect the attribute $d_{i}$ in image $x$. The model in Eq. 1 provides a
degree of explainability because a prediction must be based in natural
language attributes $d_{i}\in\mathcal{D}(c)$. At inference, we perform multi-
class classification by scoring Eq. 1 for each class, and picking the highest
scoring one, i.e. $\arg\max_{c}f_{c}(x;\mathcal{D})$.

Concept bottleneck models have traditionally been challenging to implement in
computer vision because we need to instantiate a discriminative set of
discrete attributes $\mathcal{D}(c)$, which pose challenges for gradient-based
optimization methods that are now ubiquitous in deep learning. Prior work has
relied on manual annotation of $\mathcal{D}(c)$ for each class [6], which does
not efficiently scale, or relied on other knowledge bases [5], which cannot
generalize to specialized categories.

###  3.2 Learning and Optimization

Given $C$ categories (without semantic labels), we want to learn
discriminative attributes $\mathcal{D}$ where $\mathcal{D}(c)$ are the
attributes for class $c$. We optimize the objective:

$\displaystyle\min_{\mathcal{D}}\;\mathbb{E}_{(x,y)}\left[\mathcal{L}\left(\hat{y},y\right)\right]\quad\textrm{for}\quad\hat{y}=\left[f_{1}(x;\mathcal{D}),\ldots,f_{C}(x;\mathcal{D})\right],$
(2)

where $\mathcal{L}$ is a loss function that measures the error of the
predictions $\hat{y}$ to the label $y$ of each training example. We use the
cross-entropy loss for $\mathcal{L}$.

Since attributes for each class are discrete and $f_{c}$ is not
differentiable, we optimize Eq. 2 with evolutionary search. We maintain a bank
of hypotheses $B$ for potential $\mathcal{D}(c)$, mutate them, and keep the
best-scoring hypotheses according to the objective function. Typically,
evolutionary search creates heuristics to mutate the attributes
$\mathcal{D}(c)\in B$, for example by randomly merging attributes together
(called crossover) or by randomly injecting new words from a vocabulary.
However, these heuristics are not efficient for two reasons. First, the search
space of natural language descriptions is large. Second, most heuristics do
not leverage patterns between attributes and their performance that could
drive the optimization to convergence rapidly.

We propose to replace the mutation step with a large language model and its
in-context learning capabilities. Given $k$ past hypotheses
$\\{\mathcal{D}_{t}(c),\ldots,\mathcal{D}_{t-k}(c)\\}$ for a class $c$, and
their loss, we “mutate” the next hypotheses through:

$\displaystyle\hat{\mathcal{D}}_{t+1}(c)=\textrm{LLM}\left(\mathcal{D}_{t}(c),\ldots,\mathcal{D}_{t-k}(c)\right)$
(3)

where $\mathcal{D}_{t}$ is an attribute set from iteration $t$. By using an
LLM, the mutated $\hat{\mathcal{D}}_{t+1}(c)$ benefit from natural language
priors, allowing us to efficiently search for descriptors that obey the
semantics and syntax of natural language. Secondly, the in-context learning
ability of the LLM means they will be able to find patterns in the past
hypothesis to guide the search towards attributes that are likely to minimize
the objective function. After mutation, we add $\hat{\mathcal{D}}_{t+1}(c)$ to
the bank of classifiers if it improves the objective and iterate.

We use the open-source Llama-2-70B-Instruct [67] for the large language model
$\textrm{LLM}(\cdot)$ and CLIP ViT-B/32 [1] as our vision-language model
$\phi(\cdot)$. The starting attributes are initialized randomly, with no prior
knowledge of the class name or prior information.

1

Input: Number of categories $C$, Loss function that scores attribute set
$\mathcal{L}(\mathcal{D})$ on the training set, Scalar hyper-parameters $N$,
$M$

Output: Discriminative attributes $\mathcal{D}_{\textrm{best}}$, where
$\mathcal{D}_{\textrm{best}}(c)$ is the set of attributes to recognize class
$c$

2

// Randomly initialize the classifier bank $B$

3 $B\leftarrow\\{\\}$

4 for  _$i=1,\ldots,N$_ do

5 for  _$c=1,\ldots,C$_ do

6 $\mathcal{D}_{i}(y)\leftarrow\textrm{randomly sampled attributes}$

7

8 $B\leftarrow B\;\cup\;\\{\mathcal{D}_{i}\\}$

9

// Learn a $\mathcal{D}$ that discriminates the training set $\\{(x,y)\\}$

10

11while  _not converged_ do

// Biased sampling of $M$ classifiers by expected loss
$\mathbb{E}\left[\mathcal{L}\right]$

12 $S=\\{\mathcal{D}_{1},...,\mathcal{D}_{M}\\}$ where $\mathcal{D}_{i}\sim
B\;\;\textrm{s.t.}\;\;p(\mathcal{D}_{i})\propto\mathbb{E}\left[\mathcal{L}(\mathcal{D}_{i})\right]$

// Sort by the expected loss

13 $S=(\mathcal{D}_{1},...,\mathcal{D}_{M})$ where
$\mathbb{E}\left[\mathcal{L}(\mathcal{D}_{i})\right]\geq\mathbb{E}\left[\mathcal{L}(\mathcal{D}_{j})\right]\;\;\forall_{i<j}$

// One step of evolutionary search for each category

14 for  _$c=1,\ldots,C$_ do

// Mutate the attributes with the LLM

15 $S_{c}=(\mathcal{D}_{i}(c):\mathcal{D}_{i}\in S)$

16 $\hat{\mathcal{D}}(c)\leftarrow\textrm{LLM}(S_{c})$

17

18 $S^{\prime}=(\mathcal{D}^{\prime}_{1},\ldots,\mathcal{D}^{\prime}_{M})$
where
$\forall_{k}\;\mathcal{D}^{\prime}_{i}(k)=\begin{cases}\hat{\mathcal{D}}(c)&\textrm{if
}c=k\\\ \mathcal{D}(k)&\textrm{otherwise}\end{cases}$

// Keep best attribute set

19 $B\leftarrow B\;\cup\;\arg\min_{\mathcal{D}\in
S^{\prime}}\\{\mathbb{E}\left[\mathcal{L}(\mathcal{D})\right]\\}$

20

21

22

23return $\arg\min_{\mathcal{D}\in
B}\\{\mathbb{E}\left[\mathcal{L}(\mathcal{D})\right]\\}$

24

Algorithm 1 Discriminative Attribute Set Learning

###  3.3 Evolutionary Search

Algorithm 1 shows the evolutionary search procedure to optimize the attribute
sets for image classification. We first initialize the classifier bank $B$
with random words to create $N$ initial hypotheses for $\mathcal{D}_{i}$.
During learning, we sample sets of attributes from $B$ to construct the in-
context examples for the mutation step, and repeat this process until
convergence. We bias the samples according to the loss of each attribute set
(normalized with the softmax operation). For each sampled set, we construct a
prompt per class by concatenating the attributes from class $c$ in increasing
order of performance. We prompt the LLM separately per class, which generates
the novel, mutated attributes, $\hat{\mathcal{D}}(c)$ for each class. We
evaluate the newly mutated attribute sets and add the mutated classifier with
the lowest loss to the classifier bank. We provide more implementation details
in the supplementary material.

Figure 3: Attribute evolution. We show examples of the attribute evolution for
both the pre-training and joint-training stages of learning. At the beginning,
the first generated set of attributes have little to do with the class, and by
the end of the joint-training, the learned attributes are specific to the
Greenleaf Manzanita.

###  3.4 Classifier Bank Initialization

During fine-grained classification, the optimization favors attributes that
are highly discriminative between visually similar classes, causing attributes
common to multiple classes to be discouraged. We want to encourage common
attributes at the beginning of learning so that the optimization discovers
class-specific details, instead of spurious unrelated attributes arising from
noise.

To achieve this, we use a pre-training strategy to first discover common
attributes, which serves as the initialization for joint multi-class training.
We implement the pre-training step by learning a binary classifier per class,
using the same evolutionary algorithm as before, but with an objective
function to separate one class from all others, including significantly
unrelated classes. We use the following objective for pre-training:

$\displaystyle\min_{\mathcal{D}}\;\mathbb{E}_{x_{p}}\left[f_{c}\left(x_{n};\mathcal{D}\right)\right]-\mathbb{E}_{x_{n}}\left[f_{c}\left(x_{p};\mathcal{D}\right)\right]$
(4)

where $x_{p}$ is a positive images of the class $c$, and $x_{n}$ denotes the
negative images of all the other classes. Crucially, we include classes
outside of the fine-grained dataset in the negative classes here as well. The
generated attributes across the first $200$ iterations become the
initialization of the set of attributes for the program bank initialization of
the joint training. We randomly initialize the attribute bank for binary-pre-
training with a large pool of attributes generated by an LLM about generic
visual categories. Fig. 3 shows an example of the initialization, and we
include further implementation details in the supplementary.

Figure 4: Qualitative Results. We show qualitative results for the iNaturalist
Lichen family and a KikiBouba dataset. The results illustrate two sample
images per class and the learned attributes. The learned attributes for the
Lichen hardly refer to color, as this is a common feature to all Lichen, and
instead focus on structural properties.

##  4 Experiments

###  4.1 Datasets

We validate our method by evaluating on two different datasets. First, we
evaluate on subsets of iNaturalist – a fine-grained classification dataset
with rare species that are rarely discussed online. We then evaluate on an
image dataset of novel concepts, created with new words that are even less
likely to be discussed online, as they are invented.

iNaturalist: iNaturalist [9] is a dataset for fine-grained species
classification. It contains images and annotations obtained from citizen
scientists for a large number of animal, plant, and fungus species. We
experiment with five different families and classify between five to six
species with each family. We chose families whose features varied in more
complex ways than color.

Lichen (fungi) has 6 species: elegant sunburst, golden-eye, slender orange-
bush, golden-hair, hooded sunburst, and maritime sunburst lichen. Wrasse
(fish) has 5 species: caribbean bluehead, six-bar, cortez-rainbow, moon, and
ornate wrasse. Wild rye (grass) has 5 species: squirreltail, bottlebrush,
quack grass, canada wild rye, and virginia wild rye. Manzanita (berry shrubs)
has 5 species: big-berry, pinemat, greenleaf, point-leaf, and pine-mat
manzanita. Bulrush (herb) also has has 5 species: dark-green, woolgrass,
panicled, rufous, and wood bulrush.

Kiki-Bouba: We also validate our method on completely novel concepts that do
not appear on the internet and to which language models lack familiarity. In a
surprising study, the Kiki-Bouba experiment [68] showed that people tend to
associate specific symbols with different sounds, even though the words and
the physical objects do not exist. We create a dataset of images corresponding
to non-existent concepts with generative models trained to generate images
from meaningless words [10]. Such a dataset with novel concepts makes it a
strong testbed for attribute discovery.

###  4.2 Baselines

We compare our method against several baselines. Every baseline that is not
our own has access to privileged information that we do not. Specifically, the
baselines with zero-shot attributes were generated by prompting GPT3 [69] with
the class name. The baselines that contain class names have an evident
advantage. Lastly, the gradient-based approach we constructed does not have
access to privileged information such as class names. Nonetheless, our
evolutionary method significantly outperforms this baseline for every dataset,
thus justifying our evolutionary approach. Our results show that our method
outperforms all of the baselines across all datasets, demonstrating that it
can learn attributes for specialized and undocumented visual categories.

Class Name (CLIP [1]): Our first baseline is the simple method of classifying
with the class name. For the iNaturalist dataset, which has both a common and
a scientific name per species, we report the best accuracy between using the
common name, the scientific name, and both names.

Classification by Description [5]: We implement the Classification by
Description method, proposed by [5], which generates zero-shot visual
attributes for a class by prompting GPT3 [69], and joins the class name to
each of the zero-shot attributes. Similarily to the “Class Name” baseline, for
the iNaturalist dataset, we report the best accuracy between using the common
name, the scientific name, and both names.

Zero-shot Attributes: We also use a variant of Classification by Description
where the class name is no longer appended to the attributes.

Gradient-based Approach: Instead of using an LLM to search for attributes, we
instead search for optimal input tokens to the text encoder of the VLM using
gradient descent. For a class, we find tokens that highly discriminate it from
other classes. The model optimizes for probability values over the complete
token list. However, since the explanations have to be tokens and not a
probability distribution over them, we enforce the probability distribution to
be more selective to fewer tokens as training progresses by using a
temperature parameter that decreases over the training period.

Table 1: We report the accuracy per dataset, for our method and baselines.  | iNaturalist | Kiki-Bouba
---|---|---
Method | Lichen | Wrasse | Wild Rye | Manzanita | Bulrush | KB1 | KB2
Zero-shot Attributes | 28.3 | 16.0 | 22.0 | 18.0 | 24.0 | 20.6 | 19.2
Class Name (CLIP [1]) | 23.3 | 32.0 | 32.0 | 26.0 | 26.0 | 38.7 | 38.8
Classification by Desc. [5] | 30.0 | 34.0 | 36.0 | 28.0 | 20.0 | 28.8 | 36.8
Gradient-based Approach | 23.3 | 20.0 | 40.0 | 20.0 | 20.0 | 16.7 | 55.6
Ours (1-prompt) | 31.6 | 24.0 | 44.0 | 40.0 | 22.0 | 50.3 | 47.8
Ours (10-prompt) | 48.3 | 44.0 | 58.0 | 58.0 | 42.0 | 79.2 | 59.4

Figure 5: Predictions. We show three different prediction examples. For each
example, we show our method’s prediction (first column), as well as
classification by description [5] (CBD)’s prediction (second column), and
CBD’s attributes for the ground truth class (third column). For each column,
we show the normalized probability per attribute. Below the input image, we
show the probability distributions across classes for both our method and CBD.
The results show that our learned attributes are more detailed and
discriminative of the species within the family, compared to the description
by classification (CBD) baseline. Furthermore, our method’s class probability
distributions tend to be more concentrated than CBD’s.

Varying Prompt Length: We report the results of our method with two different
prompt lengths: one and ten. In the former, the LLM prompt only has one
example of a set of attributes for a particular class, and therefore doesn’t
see the increasingly better sets of attributes. We notice a significant drop
in accuracy when in-context learning is prohibited.

Engineered Text Templates: For each method, we report the best accuracy per
dataset between the accuracy computed using scores averaged over the
engineered text templates proposed in CLIP [1], and without averaging over the
engineered text templates. We outperform all baselines with or without
engineered text templates, and report the full numbers in the supplementary.

Figure 6: Comparison of Attributes by Method. We show qualitative examples of
our learned attributes, classification by description’s attributes (CBD), and
our gradient-based approach attributes. CBD often produces reasonable
attributes, but they are not discriminative, resulting in poor recognition
accuracy. Gradient-based methods often produce poor attributes due to
optimization difficulties.

###  4.3 Fine-Grained Classification (iNaturalist)

Quantitative Result Discussion: We present the performance on accuracy for our
method and baselines in Tab. 1. Crucially, our method starts out with zero
prior knowledge on any information as to what is in a class of images, beyond
the fact that those images are grouped together. The initialized best
classifier at the beginning of optimization has no correspondence to the class
name or quality, and it is only through optimization that our method learns to
discover interpretable and reflective attributes. The first three baselines of
zero-shot descriptors, class name, and class name with zero-shot descriptors
on the other hand have prior knowledge on the class from the very start, and
therefore are at a significant advantage as they wouldn’t work for any
concepts that the internet isn’t already familiar with.

Qualitative Result Discussion: In Fig. 4, we show the fine-grained
classification results on the Lichen family. We observe that while all the
Lichen are orange and yellow in color, the only references to color within the
attributes are when the color diverges from the mean color, i.e. “reddish
tint”. The majority of the attributes reflect the Lichen’s structure, which is
due to the fact that all the Lichen are roughly the same color and grow in
similar environments, the principal discriminating feature is the structure.
In Fig. 6, we compare our learned attributes to the zero-shot attribute in
CBD, as well as the learned attributes with our gradient-based approach. The
gradient-based approach has far less descriptive and interpretable attributes.
We provide more such qualitative examples of the three methods in the
supplementary material.

In Fig. 5, we explicitly visualize the predictions across for three of the
classifiers, one per row. Each example illustrates how each method compares to
the top-performing baseline, classification by description (CBD). We show our
method’s prediction as well as CBD’s prediction, along with the relative
contribution of each attribute to the mean score. We additionally denote CBD’s
attributes for the ground truth class, along with the relative contribution to
the ground truth mean score, to visualize why CBD may have incorrectly
predicted the class. Across the three examples, we notice that CBD’s
attributes are less specific and detailed compared to ours, and that there are
many shared attributes across different classes. Shared attributes are not
useful in fine-grained classification, since the goal is to discriminate
between classes.

###  4.4 Novel Objects (Kiki-Bouba)

Quantitative Result Discussion: We outperform all baselines on the two
KikiBouba datasets. We notice that the gradient-based baseline has a large
variation in performance, and we hypothesize that this is due to the intra-
class variation changing quite significantly across classes. Across all
experiments, the gradient-based approach produces poor attributes due to
optimization difficulties, as illustrated in Tab. 1.

Qualitative Result Discussion: For the novel image results on the first
dataset of Kiki-Bouba. In Fig. 4, we show qualitative results of the learned
attribute sets for the second Kiki-Bouba dataset. We notice that compared to
the iNaturalist datasets, there is less similarity between classes, and the
discovered attributes are more object-oriented. We suspect these observations
are not unrelated, and that when an object name can be used to discriminate
between classes, the vision-language model scores highly with it.

Figure 7: Confidence. We show the mean margin of the scores between the top
and runner-up prediction, which measures the typical confidence of each model.

###  4.5 Analysis of Learning

Margin Metric: In Fig. 5, below the input image, we show the probability
distributions across classes for both our method and CBD. We observe that in
general, the probability distribution across classes for CBD is more uniform
than ours. We investigate this observation by measuring the _margin_ metric
per dataset, per method. The margin is defined as $\max_{c\in
C}f_{c}(x_{i};\mathcal{D})-\max_{\bar{c}\notin
C}f_{\bar{c}}(x_{i};\mathcal{D})$, which measures the difference in scores
between the top prediction and the runner-up, indicating the prediction
confidence [70]. We plot the margin for each method, per dataset, in Fig. 7.
The results illustrate that our method is on average twice as confident as the
other methods.

Figure 8: Dataset Bias. The squirreltail species is the only species that
commonly lives in drought habitats amongst the family, and the learned
attributes are names of plants that live in the desert. The ability to
explicitly audit bias is an advantage of our interpretable method.

Attribute Evolution: In Fig. 3, we show examples of the iterations for pre-
training, followed by the joint training. At initialization, the first sampled
set of attributes for Greenleaf Manzanita hardly relate to the class. The
closest attributes are “dark green leaves” and “green coloration”, with the
other attributes having nothing in common with images. However, by the end of
the attribute evolution, the color “green” is no longer part of the
attributes, as the other Manzanita also have green leaves, making it non-
discriminative. The final attributes contain descriptors that are particularly
descriptive of the Greenleaf Manzanita. We share more examples of attribute
evolution in the supplementary.

###  4.6 Auditing Dataset Bias

By having explicitly interpretable attributes as the bottleneck for
classification, we can directly observe whether the classifier is picking up
on dataset bias to perform prediction. An example of dataset bias can be seen
in Fig. 8. This is unique to our method, as classifiers with no concept
bottlenecks have no way of converting dataset bias into language, and previous
work in concept bottlenecks do not discover the attributes from the data
itself.

##  5 Discussion

Societal Impacts and Limitations. Explainable vision systems have the
potential to have significant practical impact, especially in specialized,
critical, and scientific domains. Such interpretable classifiers can establish
trust, as they provide insight into how a classifier reached a decision. They
allow people to audit the decision-making process, which is important for many
practical cases. Finally, they also impact education, as the classifiers can
report to a person the visual differences it has discovered, thereby helping
the person learn about the recognized concept too. However, since our approach
uses open-source LLMs, our method inherits known limitations about LLMs in
bias and inappropriate generations. As research in LLMs advance, we expect our
framework to improve too.

Conclusion. We propose a framework that integrates large language models and
evolutionary search in order to learn interpretable, discrete attributes for
visual recognition. In multiple datasets, our method outperforms existing
baselines significantly.

## References

  * [1] A. Radford, J. W. Kim, C. Hallacy, A. Ramesh, G. Goh, S. Agarwal, G. Sastry, A. Askell, P. Mishkin, J. Clark, et al., “Learning transferable visual models from natural language supervision,” in ICML, pp. 8748–8763, PMLR, 2021.
  * [2] C. Schuhmann, R. Beaumont, R. Vencu, C. Gordon, R. Wightman, M. Cherti, T. Coombes, A. Katta, C. Mullis, M. Wortsman, et al., “Laion-5b: An open large-scale dataset for training next generation image-text models,” Advances in Neural Information Processing Systems, vol. 35, pp. 25278–25294, 2022.
  * [3] V. Ferrari and A. Zisserman, “Learning visual attributes,” Advances in neural information processing systems, vol. 20, 2007.
  * [4] A. Farhadi, I. Endres, D. Hoiem, and D. Forsyth, “Describing objects by their attributes,” in 2009 IEEE conference on computer vision and pattern recognition, pp. 1778–1785, IEEE, 2009.
  * [5] S. Menon and C. Vondrick, “Visual classification via description from large language models,” ICLR, 2023.
  * [6] P. W. Koh, T. Nguyen, Y. S. Tang, S. Mussmann, E. Pierson, B. Kim, and P. Liang, “Concept bottleneck models,” in ICML, 2020.
  * [7] K. Roth, J. M. Kim, A. Koepke, O. Vinyals, C. Schmid, and Z. Akata, “Waffling around for performance: Visual classification with random words and broad concepts,” CoRR, 2023.
  * [8] J. H. Holland, Adaptation in natural and artificial systems: an introductory analysis with applications to biology, control, and artificial intelligence.  MIT press, 1992.
  * [9] G. Van Horn, E. Cole, S. Beery, K. Wilber, S. Belongie, and O. Mac Aodha, “Benchmarking representation learning for natural world image collections,” in Proceedings of the IEEE/CVF conference on computer vision and pattern recognition, pp. 12884–12893, 2021.
  * [10] M. Alper and H. Averbuch-Elor, “Kiki or bouba? sound symbolism in vision-and-language models,” Advances in Neural Information Processing Systems, vol. 36, 2024.
  * [11] S. Huang, Z. Xu, D. Tao, and Y. Zhang, “Part-stacked cnn for fine-grained visual categorization,” in CVPR, 2016.
  * [12] B. Zhou, Y. Sun, D. Bau, and A. Torralba, “Interpretable basis decomposition for visual explanation,” in ECCV, 2018.
  * [13] L. Tang, D. Wertheimer, and B. Hariharan, “Revisiting pose-normalization for fine-grained few-shot recognition,” in CVPR, 2020.
  * [14] C. H. Lampert, H. Nickisch, and S. Harmeling, “Attribute-based classification for zero-shot visual object categorization,” 2013.
  * [15] A. Frome, G. S. Corrado, J. Shlens, S. Bengio, J. Dean, M. Ranzato, and T. Mikolov, “Devise: A deep visual-semantic embedding model,” in NeurIPS, 2013.
  * [16] Z. Akata, S. Reed, D. Walter, H. Lee, and B. Schiele, “Evaluation of output embeddings for fine-grained image classification,” in CVPR, 2015.
  * [17] B. Romera-Paredes and P. Torr, “An embarrassingly simple approach to zero-shot learning,” in ICML, 2015.
  * [18] E. Kodirov, T. Xiang, and S. Gong, “Semantic autoencoder for zero-shot learning,” in CVPR, 2017.
  * [19] R. R. Selvaraju, M. Cogswell, A. Das, R. Vedantam, D. Parikh, and D. Batra, “Grad-cam: Visual explanations from deep networks via gradient-based localization,” in ICCV, 2017.
  * [20] H. Chefer, S. Gur, and L. Wolf, “Generic attention-model explainability for interpreting bi-modal and encoder-decoder transformers,” in ICCV, pp. 397–406, 2021.
  * [21] R. Fong, M. Patrick, and A. Vedaldi, “Understanding deep networks via extremal perturbations and smooth masks,” in ICCV, 2019.
  * [22] B. Zhou, A. Khosla, A. Lapedriza, A. Oliva, and A. Torralba, “Learning deep features for discriminative localization,” in CVPR, 2016.
  * [23] V. Petsiuk, A. Das, and K. Saenko, “Rise: Randomized input sampling for explanation of black-box models,” CoRR, 2018.
  * [24] V. Shitole, F. Li, M. Kahng, P. Tadepalli, and A. Fern, “One explanation is not enough: structured attention graphs for image classification,” NeurIPS, 2021.
  * [25] K. Simonyan, A. Vedaldi, and A. Zisserman, “Deep inside convolutional networks: Visualising image classification models and saliency maps,” CoRR, 2013.
  * [26] Y. Goyal, Z. Wu, J. Ernst, D. Batra, D. Parikh, and S. Lee, “Counterfactual visual explanations,” in ICML, 2019.
  * [27] V. Prabhu, S. Yenamandra, P. Chattopadhyay, and J. Hoffman, “Lance: Stress-testing visual models by generating language-guided counterfactual images,” NeurIPS, 2024.
  * [28] S. Vandenhende, D. Mahajan, F. Radenovic, and D. Ghadiyaram, “Making heads or tails: Towards semantically consistent visual counterfactuals,” in ECCV, 2022.
  * [29] P. Wang and N. Vasconcelos, “Scout: Self-aware discriminant counterfactual explanations,” in CVPR, 2020.
  * [30] C.-K. Yeh, J. Kim, I. E.-H. Yen, and P. K. Ravikumar, “Representer point selection for explaining deep neural networks,” NeurIPS, 2018.
  * [31] C.-P. Tsai, C.-K. Yeh, and P. Ravikumar, “Sample based explanations via generalized representers,” CoRR, 2023.
  * [32] Y. Sui, G. Wu, and S. Sanner, “Representer point selection via local jacobian expansion for post-hoc classifier explanation of deep neural networks and ensemble models,” in NeurIPS, 2021.
  * [33] G. Pruthi, F. Liu, M. Sundararajan, and S. Kale, “Estimating training data influence by tracking gradient descent,” CoRR, 2020.
  * [34] P. W. Koh and P. Liang, “Understanding black-box predictions via influence functions,” in ICML, 2017.
  * [35] A. Silva, R. Chopra, and M. C. Gombolay, “Cross-loss influence functions to explain deep network representations,” in AISTATS, 2020.
  * [36] H. Guo, N. Rajani, P. Hase, M. Bansal, and C. Xiong, “Fastif: Scalable influence functions for efficient model interpretation and debugging,” CoRR, 2020.
  * [37] Y. Gandelsman, A. A. Efros, and J. Steinhardt, “Interpreting clip’s image representation via text-based decomposition,” 2023.
  * [38] A. Dravid, Y. Gandelsman, A. A. Efros, and A. Shocher, “Rosetta neurons: Mining the common units in a model zoo,” in CVPR, 2023.
  * [39] M. D. Zeiler and R. Fergus, “Visualizing and understanding convolutional networks,” in ECCV, 2014.
  * [40] J. Frankle and M. Carbin, “The lottery ticket hypothesis: Training pruned neural networks,” CoRR, 2018.
  * [41] A. M. Nguyen, A. Dosovitskiy, J. Yosinski, T. Brox, and J. Clune, “Synthesizing the preferred inputs for neurons in neural networks via deep generator networks,” in NeurIPS, 2016.
  * [42] C. Jia, Y. Yang, Y. Xia, Y.-T. Chen, Z. Parekh, H. Pham, Q. Le, Y.-H. Sung, Z. Li, and T. Duerig, “Scaling up visual and vision-language representation learning with noisy text supervision,” in International conference on machine learning, pp. 4904–4916, PMLR, 2021.
  * [43] M. Cherti, R. Beaumont, R. Wightman, M. Wortsman, G. Ilharco, C. Gordon, C. Schuhmann, L. Schmidt, and J. Jitsev, “Reproducible scaling laws for contrastive language-image learning,” in CVPR, 2023.
  * [44] M. Tsimpoukelli, J. L. Menick, S. Cabi, S. Eslami, O. Vinyals, and F. Hill, “Multimodal few-shot learning with frozen language models,” NeurIPS, 2021.
  * [45] J.-B. Alayrac, J. Donahue, P. Luc, A. Miech, I. Barr, Y. Hasson, K. Lenc, A. Mensch, K. Millican, M. Reynolds, et al., “Flamingo: a visual language model for few-shot learning,” NeurIPS, 2022.
  * [46] Z. Wang, J. Yu, A. W. Yu, Z. Dai, Y. Tsvetkov, and Y. Cao, “Simvlm: Simple visual language model pretraining with weak supervision,” CoRR, 2021.
  * [47] K. Desai and J. Johnson, “Virtex: Learning visual representations from textual annotations,” in CVPR, 2021.
  * [48] J. Chen, H. Guo, K. Yi, B. Li, and M. Elhoseiny, “Visualgpt: Data-efficient adaptation of pretrained language models for image captioning,” in CVPR, 2022.
  * [49] L. H. Li, M. Yatskar, D. Yin, C.-J. Hsieh, and K.-W. Chang, “Visualbert: A simple and performant baseline for vision and language,” CoRR, 2019.
  * [50] X. Xu, C. Wu, S. Rosenman, V. Lal, W. Che, and N. Duan, “Bridgetower: Building bridges between encoders in vision-language representation learning,” in AAAI, 2023.
  * [51] H. Tan and M. Bansal, “Lxmert: Learning cross-modality encoder representations from transformers,” CoRR, 2019.
  * [52] J. Lu, D. Batra, D. Parikh, and S. Lee, “Vilbert: Pretraining task-agnostic visiolinguistic representations for vision-and-language tasks,” NeurIPS, 2019.
  * [53] Y. Zeng, X. Zhang, and H. Li, “Multi-grained vision language pre-training: Aligning texts with visual concepts,” CoRR, 2021.
  * [54] L. Yuan, D. Chen, Y.-L. Chen, N. Codella, X. Dai, J. Gao, H. Hu, X. Huang, B. Li, C. Li, et al., “Florence: A new foundation model for computer vision,” CoRR, 2021.
  * [55] L. H. Li, P. Zhang, H. Zhang, J. Yang, C. Li, Y. Zhong, L. Wang, L. Yuan, L. Zhang, J.-N. Hwang, et al., “Grounded language-image pre-training,” in CVPR, 2022.
  * [56] L. Li, Z.-Y. Dou, N. Peng, and K.-W. Chang, “Desco: Learning object recognition with rich language descriptions,” NeurIPS, 2024.
  * [57] J. Li, D. Li, C. Xiong, and S. Hoi, “Blip: Bootstrapping language-image pre-training for unified vision-language understanding and generation,” in ICML, 2022.
  * [58] J. Li, D. Li, S. Savarese, and S. Hoi, “Blip-2: Bootstrapping language-image pre-training with frozen image encoders and large language models,” CoRR, 2023.
  * [59] R. Mokady, A. Hertz, and A. H. Bermano, “Clipcap: Clip prefix for image captioning,” CoRR, 2021.
  * [60] Z. Luo, Y. Xi, R. Zhang, and J. Ma, “A frustratingly simple approach for end-to-end image captioning,” CoRR, 2022.
  * [61] S. Pratt, I. Covert, R. Liu, and A. Farhadi, “What does a platypus look like? generating customized prompts for zero-shot image classification,” in Proceedings of the IEEE/CVF International Conference on Computer Vision, pp. 15691–15701, 2023.
  * [62] A. Yan, Y. Wang, Y. Zhong, C. Dong, Z. He, Y. Lu, W. Y. Wang, J. Shang, and J. McAuley, “Learning concise and descriptive attributes for visual recognition,” in CVPR, 2023.
  * [63] B. Romera-Paredes, M. Barekatain, A. Novikov, M. Balog, M. P. Kumar, E. Dupont, F. J. Ruiz, J. S. Ellenberg, P. Wang, O. Fawzi, et al., “Mathematical discoveries from program search with large language models,” Nature, 2023.
  * [64] S. Yu, S. Liu, Z. Lin, D. Pathak, and D. Ramanan, “Language models as black-box optimizers for vision-language models,” CoRR, 2023.
  * [65] S. Jin, X. Jiang, J. Huang, L. Lu, and S. Lu, “Llms meet vlms: Boost open vocabulary object detection with fine-grained descriptors,” CoRR, 2024.
  * [66] S. Han, L. Zhuo, Y. Liao, and S. Liu, “Llms as visual explainers: Advancing image classification with evolving visual descriptions,” CoRR, 2023.
  * [67] H. Touvron, L. Martin, K. Stone, P. Albert, A. Almahairi, Y. Babaei, N. Bashlykov, S. Batra, P. Bhargava, S. Bhosale, et al., “Llama 2: Open foundation and fine-tuned chat models,” CoRR, 2023.
  * [68] V. S. Ramachandran and E. M. Hubbard, “Synaesthesia–a window into perception, thought and language,” Journal of consciousness studies, vol. 8, no. 12, pp. 3–34, 2001.
  * [69] B. Mann, N. Ryder, M. Subbiah, J. Kaplan, P. Dhariwal, A. Neelakantan, P. Shyam, G. Sastry, A. Askell, S. Agarwal, et al., “Language models are few-shot learners,” CoRR, 2020.
  * [70] P. Melville and R. J. Mooney, “Diverse ensembles for active learning,” in Proceedings of the twenty-first international conference on Machine learning, p. 74, 2004.


# DMseg: a Python algorithm for de novo detection of differentially or
variably methylated regions

Xiaoyu Wang  Fred Hutchinson Cancer Center, Public Health Sciences Division,
Seattle, 98109, USA  Ming Yu  Fred Hutchinson Cancer Center, Clinical Research
Division, Seattle, 98109, USA  William Grady  Fred Hutchinson Cancer Center,
Public Health Sciences Division, Seattle, 98109, USA  Fred Hutchinson Cancer
Center, Clinical Research Division, Seattle, 98109, USA  University of
Washington, Department of Medicine, School of Medicine, Seattle, 98195, USA
Ziding Feng  Fred Hutchinson Cancer Center, Public Health Sciences Division,
Seattle, 98109, USA  University of Washington, Department of Biostatistics,
School of Public Health, Seattle, 98195, USA  Wei Sun  Fred Hutchinson Cancer
Center, Public Health Sciences Division, Seattle, 98109, USA  University of
North Carolina, Department of Biostatistics, Chapel Hill, 27516, USA
University of Washington, Department of Biostatistics, School of Public
Health, Seattle, 98195, USA  James Y Dai  Fred Hutchinson Cancer Center,
Public Health Sciences Division, Seattle, 98109, USA  University of
Washington, Department of Biostatistics, School of Public Health, Seattle,
98195, USA  Correspondence<EMAIL_ADDRESS>
###### Abstract

Detecting and assessing statistical significance of differentially methylated
regions (DMRs) is a fundamental task in methylome association studies. While
the average differential methylation in different phenotype groups has been
the inferential focus, methylation changes in chromosomal regions may also
present as differential variability, i.e., variably methylated regions (VMRs).
Testing statistical significance of regional differential methylation is a
challenging problem, and existing algorithms do not provide accurate type I
error control for genome-wide DMR or VMR analysis. No algorithm has been
publicly available for detecting VMRs. We propose DMseg, a Python algorithm
with efficient DMR/VMR detection and significance assessment for array-based
methylome data, and compare its performance to Bumphunter, a popular existing
algorithm. Operationally, DMseg searches for DMRs or VMRs within CpG clusters
that are adaptively determined by both gap distance and correlation between
contiguous CpG sites in a microarray. Levene test was implemented for
assessing differential variability of individual CpGs. A likelihood ratio
statistic is proposed to test for a constant difference within CpGs in a DMR
or VMR to summarize the evidence of regional difference. Using a stratified
permutation scheme and pooling null distributions of LRTs from clusters with
similar numbers of CpGs, DMseg provides accurate control of the type I error
rate. In simulation experiments, DMseg shows superior power than Bumphunter to
detect DMRs. Application to methylome data of Barrett’s esophagus and
esophageal adenocarcinoma reveals a number of DMRs and VMRs of biological
interest.

DMseg is available at https://pypi.org/project/DMseg

## Introduction

Epigenetics, DNA modifications that do not change the underlying sequence, may
provide an interface between environmental insult, genetic susceptibility and
disease development. There is increasing interest to study epigenetics and
understand the etiology of complex diseases such as cancers [1, 2]. The most
common epigenetic markers are CpG methylation, a chemical modification of a
cytosine (methyl-cytosine) that is immediately followed by a guanine. With
advances in high-throughput technologies, it is now possible to measure DNA
methylation for nearly all 28 million CpGs in the human genome by whole-genome
bisulfite squencing (WGBS), or more economically, a highly selected subset
with direct functional relevance measured by microarrays such as Illumina
Infinium HumanMethylation450 and EPIC BeadChips. Methylome association studies
are now routinely conducted to understand etiology or identify disease
biomarkers [3].

In human genome, CpG methylation sites are irregularly spaced, with blocks of
correlated CpGs, e.g., promoters, CpG island and CpG island shores.
Consequently, when comparing samples under different biological conditions,
differential methylation sites were often found in regions with contiguous
CpGs simultanuously hyper- or hypo- methylated [4]. The identification and
characterisation of differentially methylated regions (DMRs) between phenotype
groups have been of prime interest in methylome association studies [5]. For
cancer methylome, differential methylation can often present as mean
methylation shift as well as increased stochastic variance, the latter of
which is due to the heterogeneous nature of cancers [6, 7]. The main
methodological challenge is to discover differentially methylated region (DMR)
or variably methylated region (VMR) in an efficient manner, and to provide an
accurate assessment of statistical significance of these regional findings.

There are a number of DMR calling algorithms for the restricted CpG coverage
in HM450 arrays and EPIC BeadChips [5, 8, 9]. Nearly all existing DMR calling
methods start from computing individual CpG associations and identifying
contiguous ones with evidence of differential methylation. Local smoothing of
probe-level statistics with respect to chromosomal coordinates is commonly
adopted, for example, Bumphunter and DMRcate. Except for Bumphunter, rigorous
evaluation of statistical significance for DMR findings has not always
conducted due to computational burden required by a permutation or bootstrap
procedure. Indeed, in our data applications it may take several hours for
Bumphunter to compute family-wise error rate for discovered DMRs. On the other
hand, methods for detecting VMR have been proposed [10], though computational
software for genome-wide VMR are not yet readily available.

In this article, we develop a unified method for detecting and assessing de
novo DMR/VMR in methylome association studies. A computationally efficient
Python algorithm, DMseg, is released at PyPi.org to implement the proposed
method. Comparing to existing methods and algorithms, DMseg made several
contributions. First, a likelihood ratio test (LRT) statistic assuming a
common difference across contiguous CpGs is used to summarize the evidence of
DMR/VMR. This test provides some level of smoothing by averaging CpG-level
associations, though much simpler than local smoothing based on chromosomal
coordinates. Second, we show that test statistics for DMRs are typically not
exchageable between DMRs with different numbers of CpGs. To unbiasely assess
the false positive error rate, a stratified permutation test was implemented
in DMseg that accounts for non-exchangeability of LRT statistic. We conduct
extensive simulations to evaluate and compare DMseg to bumphunter, the only
existing algorithm with significance assessment for DMR.

##  Methods & Datasets

The objective of DMseg is to enable fast detection of DMR and VMR for
microarray-based methylome association studies, and to provide an accurate
assessment of statistical significance of the resultant DMRs and VMRs. Many
algorithms exist for detecting DMR [5, 8, 9], though few is capable of
conducting a rigorous hypothesis test for DMR findings.

We use methylome studies on progression of a cancer precursor lesion,
Barrett’s esophagus, to esophageal adenocarcinoma to illustrate the motivation
and the utility of DMseg. The scientific goal is to discover and validate
early detection methylation biomarkers for malignant transformation of BE, so
that preventative measures can be undertaken to intercept cancer progression.
Two methylome datasets were downloaded from GEO: the first dataset contains 64
normal squamous esophagus samples, 19 Barrett$\prime$s esophagus (BE) and 125
esophageal adenocarcinoma (EAC) samples from Australia [11]. After quality
control, there were 372376 CpGs for analysis; the second set is from 33 normal
squamous samples, 59 BE samples from cancer free patients, and 23 EAC samples
from United States [12]. All samples were fresh frozen and profiled by
Infinium Human Methylation 450K BeadChips.

### Algorithm modules and workflow

The input datasets for DMseg include an epidemiologic dataset containing the
primary comparison groups for association and adjusting covariates, and a
methylation data matrix for processed methylome data. Suppose there are $n$
samples with DNA methylome data available for the association analysis,
denoted by $(\mbox{\bf Y}_{i},X_{i},W_{i})$, for $i=1,...,n$, where $\mbox{\bf
Y}_{i}$ is a length $p$ vector for methylation beta values or M-values in $p$
CpG sites, $X_{i}$ is the group indicator of the primary interest, for example
the cancer ($X_{i}=1$) or normal sample ($X_{i}=0$), and $W_{i}$ is the vector
of additional covariates such as age and gender that should be adjusted for.
Let $Y_{ij}$ denote the beta value or M-value of the $j^{th}$ CpG site for the
$i^{th}$ sample.

Algorithmically, DMseg is composed of four modules (Figure 11): determine
clusters of CpGs by distance and correlation; compute CpG-level association
statistics; search clusters for candidate DMRs and computing likelihood ratio
statistics; assess significance of DMRs by a permutation or bootstrap
procedure.

Figure 1: The modules and the workflow of DMseg algorithm

Determining clusters of CpGs:  The CpG sites selected for the HM450 array and
the EPIC BeadChip represent a small subset (1.5-3%) of the entire 28 million
CpG sites in the human genome, mostly distributed in the CpG-rich regions such
as promoters, CpG islands and shores. The spacing between CpGs are highly
irregular: clusters of densely distributed, highly correlated CpGs are
separated by long intergenic regions with scarce CpGs. Allowing parallel
searching, a computationally efficient approach is to first identify clusters
and then detect DMRs within clusters, as implemented in bumphunter. DMseg
determines the two contiguous CpGs belong to a cluster if the distance between
them is less than a user defined number of base pairs (e.g. 500), or the
correlation between contiguous CpGs greater than 0.6 (the median correlation
coefficent among contiguous CpGs within the clusters defined by gap distance
only). This criterion is different from bumphunter [5], which only uses the
gap distance to determine clusters. The motivation is that correlation may not
always decay linearly with chromosomal distance in base pairs, and there are
CpGs far away in chromosomal coordinates yet highly correlated.

Using the HM450 data for BE/EAC samples as an example, we obtained 170156
clusters (including 1-CpG clusters) by the maximum gap distance 500 bp. Adding
the correlation criterion will reduce the cluster number to 150676, of which
19480 clusters were connected due to correlation. Removing 1-CpG clusters will
result in 54688 clusters, from which DMRs were searched for.

Computing CpG-level association statistics:  Users can choose whether to use
methylation beta values or M values for calculating individiual CpG
association for testing differentially methylated CpGs. The former yields
regression coefficents are directly interpretable as the average percentage of
methylation change between comparison groups, while the latter works better
for statistical performance of linear models. Let $Y_{ij}$ denote the
methylation value for the $i^{th}$ sample and the $j^{th}$ CpG site. For
testing differentially methylated CpGs, the following linear model is fitted
for each of CpGs in the array,

$\displaystyle\mathbb{E}(Y_{ij})$ $\displaystyle=$
$\displaystyle\beta_{0j}+\beta_{1j}X_{i}+\beta_{2j}W_{i},$ (1)

where $\beta_{1j}$ is the association parameter of interest. For testing
variably methylated CpGs, $Y_{ij}$ is the methylation M-value and the
following regression model will be fitted

$\displaystyle\mathbb{E}(|Y_{ij}-\widetilde{m}_{g(i)j}|)$ $\displaystyle=$
$\displaystyle\alpha_{0j}+\alpha_{1j}X_{i}+\alpha_{2j}W_{i},$ (2)

where $g(i)$ is the group (e.g., tumor or normal) label for $i^{th}$ sample,
$\widetilde{m}_{g(i)j}$ is the sample median for $j^{th}$ CpG site in $g(i)$
group, $|Y_{ij}-\widetilde{m}_{g(i)j}|$ is the absolute difference to the
corresponding group median. Model (2) implements the Levene test for
differential variability [13]. Because the predictors and the design matrix
for each of CpG association models are identical, the estimated $\beta_{1j}$s
and its estimated standard errors can be computed by sharing common
intermediary elements of ordinary least squares regression for different CpGs.
Therefore the algorithm does not need to iterate through CpGs one by one,
which substantially reduces the computation time.

Detecting DMR/VMR and computing likelihood ratio statistics:  We consider
candidate DMR/VMR to be regions with at least $k$ contiguous CpGs that are
differentially methylated, where $k=2$ or 3 as defined by users.
Operationally, searching for candidate DMR/VMR is conducted within clusters
that are previously determined, using $z$-score for individual CpG association
$\geq$ 1.96 as a soft threshold for differential methylation. Two features of
DMseg permit flexibility in accommodating diverse genomic irregularities:
first, DMseg allows switching sign of CpG association (hypermethylation or
hypomethylation) within a DMR, similar to DMRcate [9], since sometimes
hypermethylation in gene bodies may also be associated with upregulation of
gene expressions; second, though DMseg don’t employ any smoothing method for
CpG level association statistics, if two adjacent DMRs or two adjacent VMRs
within one cluster are separated by merely a single CpG, and that CpG shows a
moderate level of association, e.g., $z$-score $\geq$ 1.64 (user can define
this softer threshold), then the two DMRs or VMRs will be merged through the
CpG.

After a candidate DMR was detected, a likelihood ratio statistic will be
computed for the DMR, comparing the null model that there is no association
for any CpGs within the DMR, to the alternative model that all CpG-level
association parameters within a DMR are a same scalar. Suppose
$\widehat{\mbox{\boldmath$\beta$}}$ =$(\hat{\beta}_{11},...,\hat{\beta}_{1L})$
are estimated association parameters for $L$ CpGs within a candidate DMR. Let
$\hat{\sigma}^{2}_{11},...,\hat{\sigma}^{2}_{1L}$ denote the variance of
$\hat{\beta}_{11},...,\hat{\beta}_{1L}$ respectively. When sample size $n$ is
sufficiently large, the distribution of
$\hat{\sigma}^{2}_{11},...,\hat{\sigma}^{2}_{1L}$ can be approximated by a
Gaussian distribution with its corresponding variance. Under the alternative
hypothesis that all $\beta_{1j}$ are the same constant $\bar{\beta}_{1}$, the
maximum likelihood estimate of $\bar{\beta}_{1}$ is $\mbox{\bf
J}^{\prime}\Sigma^{-1}\widehat{\mbox{\boldmath$\beta$}}$/$\mbox{\bf
J}^{\prime}\Sigma^{-1}\mbox{\bf J}$, where J is the length-$L$ vector of 1,
and $\Sigma$ is $L\times L$ diagonal matrix with its diagonal elements
$\hat{\sigma}^{2}_{11},...,\hat{\sigma}^{2}_{1L}$. The LRT statistic is
therefore

$\widehat{\mbox{\boldmath$\beta$}}^{\prime}\Sigma^{-1}\widehat{\mbox{\boldmath$\beta$}}-(\widehat{\mbox{\boldmath$\beta$}}-\bar{\beta}_{1})^{\prime}\Sigma^{-1}(\widehat{\mbox{\boldmath$\beta$}}-\bar{\beta}_{1}).$
(3)

The larger this LRT is, the greater evidence for a significant DMR. As we show
next, a larger LRT could be resulted from a larger difference in mean
methylations, or a DMR with more CpGs included.

Figure 22 shows an example DMR in chromosome 8. There are 22 CpGs within this
cluster (starting position 145728203, ending position 145729799), all of which
showed hypermethylation in esophageal adenocarcinoma samples relative to
Barrett’s esophagus samples. The z-scores of all 22 CpG associations are
greater than 2.0, suggesting a highly significant DMR. DMseg fitted a
horizontal line to the 22 differences of beta values with its height slightly
greater than 0.1, and the LRT for this DMR is 599.5.

Figure 2: An example of DMR fitted by DMseg algorithm

Significance test of DMRs by permutation or bootstrap:  One distinguishing
feature of DMseg is to conduct a permutation/bootstrap test for assessing the
significance of DMR. This can be computationally demanding for genome-wide
testing. For example, it takes several hours for bumphunter to conduct its
bootstrap-based testing procedure. Moreover, hypothesis testing for a data-
driven DMR finding can be challenging because: i) it is difficult to establish
the null distribution of such data-driven DMR; ii) the test statistic, e.g.,
LRT for DMseg may not be “exchangeable” for different clusters. Figure 33
presents the permutation null distributions from three clusters with different
sizes (5 CpGs, 22 CpGs, 97 CpGs). Clearly clusters with a greater number of
CpGs have larger LRT values: 94.5% of permutations from the 5-CpG cluster
yield no DMR finding (therefore LRT=0), and the 95% quantile of its null
distribution (including LRT=0 permutations) is 8.6; 87% of permutations from
the 22-CpG cluster yield no DMR finding, and the 95% quantile of its null
distribution is 20.4; 59.3% of permutations from the 97-CpG cluster yield no
DMR finding, and the 95% quantile of its null distribution is 59.6. In
summary, more CpGs in the cluster to start the searching lead to a higher
probability of finding DMR, and a DMR with more CpGs yields a larger LRT.

Figure 3: Clusters with different number of CpGs have different null
distributions for LRT.

This “non-exchageability” property for LRT statistics suggests that one cannot
pool null LRT statistics across clusters with different numbers of CpGs.
Ideally, a large number of permutations need to be conducted for each cluster
to result in a p-value with $10^{-8}$ resolution, which is computationally
prohibitive. DMseg implements a stratified pooling strategy: pooling
permutation LRTs as the null distribution for computing p-values from 4 strata
with CpG numbers in (0,10],(10,20],(30,40] and more than 40. For clusters
without DMR findings, their p-values are set to be 1. FWERs for top DMR
candidates are then calculated from these p-values.

## Results

### Synthetic data to evaluate type I error rate and power

A simulation experiment was conducted to evaluate the type I error rate and
the power of DMseg and Bumphuner. To mimic real methylome data, the synthetic
data were created using on the normal DNA methylation data in the BE/EAC
example [11] as the backbone, and adding a gradient of effect sizes for
differential methylation. Specifically, the data for the null hypothesis were
constructed by randomly splitting the 64 normal squamous esophagus samples to
two groups for 100 times, each with 32 samples. To reduce the computation
time, we restricted to all CpGs data in chromosome 10, and supplemented CpGs
from clusters with numbers of CpGs $>$ 20 from all other chromosomes. The
rationale is to test the performance in the full spectrum of CpG clusters,
particularly large clusters with many CpGs. In total, there are 29165 CpGs
included, and 9975 CpGs of them are from other chromosomes. These CpGs form
3197 clusters.

To generate a DMR and test the power, we selected 10 CpGs in the _Vimentin_
gene and added signals in either the beta value scale or the M value scale to
every sample in the designated “case” group. _Vimentin_ gene in Chromosome 10
is a known methylation biomarker for detecting BE and EAC [14]. Figure 44
shows the mean beta values among 10 CpGs (from cg05151811 to cg02746869) in
_VIM_ gene for BE, EAC, and normal samples. There are marked differences in
mean beta values between normal samples and BE/EAC samples for the 10 CpGs in
the promoter region, though the differences between BE and EAC are less
pronounced. As shown in Table 1, to evaluate the full-range power performance
four levels of DMR effect sizes were added to the beta values of the 10 CpGs
for every sample in the “case” group, and similarly five levels of effect
sizes in the M-value scale. Note that when effect size is set to be zero, the
family wise error rate (FWER, the probability of declaring one or more false
positives) was evaluated.

The comparator Bumphunter can be used with customized options. For example,
users can pick a cutoff value for the minimal differences of mean beta values
or mean M values for contiguous CpGs to be considered as DMR, or users can let
Bumphunter decide the cutoff value by simulations. While Bumphunter proposed
to use smoothed CpG methylation differences for detecting DMR, users can
choose the option of no smoothing and therefore using the original methylation
differences. The results for Bumphunter under these different options were
shown in Table 1. In either scale, DMseg preserved the correct FWER in 100
permutation datasets, and delivered a much improved power over Bumphunter,
regardless of the options used in Bumphunter. When Bumphunter was set to
automatically select a cut-off value for detecting DMR, FWER is inflated
whether smoothing was used or not. The smoothing option clearly worsened the
performance of Bumphunter

Table 2 shows the results when only the last five CpGs were used in the
simulation experiment. This scenario was created to evaluate power in weaker
signal-to-noise settings, and may present challenges to the smoothing option
of Bumphunter. The superior performance of DMseg over Bumphunter is more
drastic for this scenario with 5 differentially methylated CpGs.

Figure 4: The mean CpG methylation levels for the ten CpGs in VIM gene chosen
for the simulation study. Table 1: Simulation results for type I error rate
and power when adding differential methylation to 10 CpGs in VIM gene

| Analysis
---
scale
| Effect size in
---
chosen scale
| Bumphunter No smoothing
---
(smoothing)
DMseg
|  | Set thresholda |  | Automatically
---
pick threshold
Beta | FWER | 0 | 0.03 (0.03) | 0.08 (0.11) | 0.05
| Power | 0.05 | 0 (0) | 0 (0) | 0.53
|  | 0.1 | 0.07 (0) | 0 (0) | 0.97
|  | 0.15 | 0.86 (0.72) | 0.05 (0) | 1
M-value | FWER | 0 | 0.05 (0.06) | 0.06 (0.07) | 0.06
| Power | 1 | 0 (0) | 0.07 (0) | 0.48
|  | 1.2 | 0 (0) | 0.29 (0) | 0.74
|  | 1.4 | 0 (0) | 0.54 (0) | 0.94
|  | 1.6 | 0 (0) | 0.72 (0.01) | 0.94

a 0.1 used for beta values, 0.28 for M-value

Table 2: Simulation results for type I error rate and power when adding
differential methylation to 5 CpGs in VIM gene

| Analysis
---
scale
| Effect size in
---
chosen scale
| Bumphunter No smoothing
---
(smoothing)
DMseg
|  | Set thresholda |  | Automatically
---
pick threshold
Beta | FWER | 0 | 0.03 (0.06) | 0.08 (0.11) | 0.05
| Power | 0.1 | 0 (0) | 0 (0) | 0.85
|  | 0.15 | 0.07 (0) | 0 (0) | 0.98
|  | 0.2 | 0.01 (0.02) | 0 (0) | 1
M-value | FWER | 0 | 0.05 (0.06) | 0.06 (0.08) | 0.06
| Power | 1.6 | 0 (0) | 0 (0) | 0.44
|  | 2 | 0 (0) | 0.06 (0) | 0.75
|  | 2.4 | 0 (0) | 0.33 (0) | 0.95
|  | 2.8 | 0 (0) | 0.55 (0.01) | 0.97

a 0.1 used for beta values, 0.28 for M-value

### HM450 data example: methylome comparison between BE and EAC

DMseg was used in a biomarker study to identify DMR and VMR between BE and
EAC. The Australian dataset was used as the discovery set [11], and the US
dataset was used as the validation set [12]. As a comparison, Bumphunter did
not yield any significant DMR at FWER 0.05 level.

Using DMseg, we identified 47 significant DMRs in 16 chromosomes and 2
significant VMRs in 2 chromosomes after correcting for multiple testing
(FWER<0.05). A linux machine with Intel(R) Xeon(R) Gold 6254 CPU (3.1 GHz),
and 64 GB memory was used in this study. The most recent version of
Bumphunter) (V1.38.0) was run on R (V4.2.0), and DMseg was run on Python
(V3.8.12). It took DMseg 12 minutes 10 seconds to run DMR analysis for this
discovery set (23 minutes 33 seconds if using Bumphunter). The top 10 DMRs and
the two VMRs are shown in Table 3, ordered by FWER. The full list of DMRs are
included in Supplementary Table 1. For each local cluster that contains a DMR
or VMR identified in the discovery set, we ran DMseg in the validation set see
if the same or similar DMR can be detected in the cluster. Sometimes not all
CpGs in a DMR were available for validation, due to removal of CpGs during the
quality control steps of the two datasets. Strikingly, all DMRs and VMRs
identified in the discovery set were validated with a p-value less than 0.05
(Table 3, Supplementary Table 1). A high degree of concordance was observed
for LRTs and segment means for DMRs in the discovery set and in the validation
set, supporting the robustness of these DMR findings.

Figure 55 shows a significant DMR in Chromosome 14 detected by DMseg, ranked
the 39th among the 47 DMRs. This is an example in which merging long separated
CpGs by correlation led to new discoveries. This is a region in _MTA1_ gene,
which have been reported to promote tumorigenesis and development of
esophageal squamous cell carcinoma, also a potential indicator for assessing
the malignant potential of colorectal and gastric carcinomas [15, 16]. The
correlation matrix of CpGs in this region is shown by a heatmap in the left
panel; The difference of mean beta along chromosome coordinates are shown in
the right panels (CpGs in different clusters are indicated by different
colors). The Cluster B is an example of connecting contiguous CpGs with high
correlation but a longer distance than the max gap 500 bp: without merging by
correlation, there would be 3 clusters in Cluster B with the sizes of 3 CpGs,
1 CpG, 3 CpGs respectively. Connecting these CpGs created a larger cluster,
which led to a statistically significant DMR.

Figure 5: A DMR detected by DMseg in the Discovery set. (A) Pearson
correlation among CpGs in Clusters A-C. Dashed boxes in different colors show
correlation among CpGs in different clusters, and the solid purple box
indicates correlation among CpGs within the DMR. (B) DMR (the purple box)
detected in this region. (C) The same DMR in the Validation set.

The two significant VMRs in Table 3 are not overlapped with any of the 47
DMRs, though they also have a moderate level of mean differences between BE
and EAC samples (p-values for testing DMR are $2.52\times 10^{-5}$ and
$7.37\times 10^{-4}$, with corresponding FWER 0.496 and 1). For example, the
top ranked VMR (Figure 66) contains 6 CpGs in Chromosome 20 in the promoter
region of _DIDO1_ gene, which has been reported to play an important role in
inhibiting gastric cancer progression [17]. All six CpGs showed low
methylation in the BE samples, and significantly increased variability and
hypermethylation (though did not reach 0.05 in FWER) in EAC samples. The mean
standard deviation of beta values for 6 CpGs in BE samples is 0.046 (mean
methylation level 0.09), and it increases to 0.25 for EAC samples (mean
methylation level 0.31). The statistical significance for the increased
variability is visibly stronger than that for hypermethylation.

Table 3: Differentially methylated regions in the BE-EAC comparison | chr | gene | Discovery set | Validation set
---|---|---|---|---
|  |  | start CpG | end CpG | #CpG | mean | LRT | p-value | FWER | #CpG | mean | LRT | p-value
DMR | chr2 | _CAPN10_ | cg08589214 | cg16615776 | 5 | 0.22 | 599.5 | <$1\times 10^{-6}$ | <0.002 | 5 | 0.21 | 326.4 | <0.002
| chr11 | _AP2A2;MUC6_ | cg21625737 | cg20839554 | 6 | 0.19 | 409.4 | <$1\times 10^{-6}$ | <0.002 | 6 | 0.20 | 275.5 | <0.002
| chr3 | _EPHB3_ | cg03867377 | cg10769844 | 5 | 0.15 | 352.5 | <$1\times 10^{-6}$ | <0.002 | 5 | 0.14 | 124.7 | <0.002
| chr11 | _NUMA1_ | cg02151120 | cg05874355 | 6 | 0.22 | 345.0 | <$1\times 10^{-6}$ | <0.002 | 7 | -0.05;0.22 | 329.2 | <0.002
| chr1 | _SNORD103A/PUM1_ | cg24722354 | cg00581541 | 4 | 0.22 | 311.7 | <$1\times 10^{-6}$ | <0.002 | 3 | 0.17 | 66.6 | <0.002
| chr13 | _A2LD1_ | cg02425595 | cg25710107 | 9 | 0.14 | 303.6 | <$1\times 10^{-6}$ | <0.002 | 9 | 0.11 | 155.2 | <0.002
| chr4 | _LOC100130872/SPON2_ | cg20756245 | cg17232217 | 16 | 0.13 | 293.9 | <$1\times 10^{-6}$ | <0.002 | 17 | 0.16 | 507.4 | <0.002
| chr1 | _GJC2;GUK1_ | cg22844623 | cg20848130 | 9 | 0.16 | 276.0 | <$1\times 10^{-6}$ | <0.002 | 9 | 0.17 | 258.4 | <0.002
| chr14 | _ZFYVE21_ | cg13154413 | cg25580656 | 6 | 0.24 | 252.1 | <$1\times 10^{-6}$ | <0.002 | 6 | 0.17 | 169.5 | <0.002
| chr14 | _JDP2_ | cg25102206 | cg22143352 | 4 | 0.24 | 247.5 | <$1\times 10^{-6}$ | <0.002 | 4 | 0.23 | 322.8 | <0.002
VMR | chr20 | _DIDO1_ | cg02461114 | cg04226366 | 6 | 1.38 | 179.4 | $4.39\times 10^{-8}$ | 0.002 | 6 | 1.48 | 432.8 | <0.002
| chr19 | _ZNF154_ | cg03142586 | cg27324426 | 10 | -0.72 | 164.3 | $1.32\times 10^{-7}$ | 0.006 | 4 | -0.43 | 20.5 | 0.032

Figure 6: A VMR detected by DMseg. (A) VMR in the Discovery set. (B) The same
VMR in the Validation set.

### Epic array data application

Other than the HM450 array, DMseg can analyze data generated by the Infinium
Methylation EPIC microarray, which provides a higher throughput platform to
quantify methylation at over 850,000 CpG sites on the human genome. To test
speed, we analyzed a Fred Hutch dataset containing 11 Gastroesophageal reflux
disease (GERD) samples and 10 BE samples. This data include 760241 CpGs after
quality control steps. Total number of clusters (including 1 CpG clusters) is
388459 before merging, 367483 after merging. These numbers are about double of
those for BE and EAC data on the HM450 platform. It took 16 minutes 41 seconds
for DMseg to run DMR analysis (34 minutes 4 seconds for Bumphunter). Due to
the small sample size, no DMR was found with FWER $<$ 0.05. On the other hand,
one VMR with 7 CpGs (from cg05087455 to cg19240637 in Chromosome 2) was found
in _RNF144A_ gene with LRT = 97.2 and FWER = 0.046.

## Discussion

In this work we proposed a fast and accurate algorithm for detecting
statistically significant DMR/VMR for array-based methylome data, DMseg. The
superior performance of DMseg can be attributed to several features of DMseg:
1) DMseg uses Z-score statistics to detect contiguous DMR/VMR, which is
statistically principled and more sensitive than the differential methylation
regression coefficients used in Bumphunter; 2) DMseg uses a likelihood ratio
test statistic that fits a constant mean to a DMR, averaging differential
methylation effect sizes among CpGs in a DMR. This avoids commonly used
smoothing procedures (loess or kernel smoothing) that requires bandwidth
selection. The latter is quite challenging for irregularity and sparse
distribution of CpG coverage in a methylation array; 3) DMseg defines clusters
within which searching for DMR is conducted using both gap distance and
correlation, which is more data-adaptive than Bumphunter. This creates longer
clusters and yields more DMRs; 4) DMseg pools null distributions of LRT for
clusters with similar numbers of CpGs, effectively increasing the power for
detecting shorter DMRs.

The most computationally expensive step of DMseg is the permutation test,
which entails calculation of differential methylation regression for genome-
wide CpGs. Long clusters may require as many as 5000 permutations to get a
high-resolution p-value. This is necessary for accurate assessment of the type
I error rate for candidate DMR findings, the distictive feature of DMseg.
Extending this algorithm to whole-genome bisulphite sequencing data and
accounting for variability of read depth is feasible, as have been done by
DMRcate [18].

## References

  * [1] Petronis, A.  Epigenetics as a unifying principle in the aetiology of complex traits and diseases.  _Nature_ 465, 721–727, DOI: 10.1038/nature09230 (2010).
  * [2] Esteller, M.  Epigenetics in cancer.  _N. Engl. J. Med_ 358, 1148–1159 (2008).
  * [3] Relton, C. L. & Davey Smith, G.  Epigenetic epidemiology of common complex disease: prospects for prediction, prevention, and treatment.  _PLoS medicine_ 7, e1000356, DOI: 10.1371/journal.pmed.1000356 (2010).
  * [4] Irizarry, R. A. _et al._ The human colon cancer methylome shows similar hypo- and hypermethylation at conserved tissue-specific cpg island shores.  _Nature genetics_ 41, 178–186, DOI: 10.1038/ng.298 (2009).
  * [5] Jaffe, A. E. _et al._ Significance analysis and statistical dissection of variably methylated regions.  _Biostatistics (Oxford, England)_ 13, 166–178, DOI: 10.1093/biostatistics/kxr013 (2012).
  * [6] Hansen, K. D. _et al._ Increased methylation variation in epigenetic domains across cancer types.  _Nature genetics_ 43, 768–775, DOI: 10.1038/ng.865 (2011).
  * [7] Phipson, B. & Oshlack, A.  Diffvar: a new method for detecting differential variability with application to methylation in cancer and aging.  _Genome biology_ 15, 465, DOI: 10.1186/s13059-014-0465-4 (2014).
  * [8] Pedersen, B. S. _et al._ Comb-p: software for combining, analyzing, grouping and correcting spatially correlated p-values.  _Bioinformatics (Oxford, England)_ 28, 2986–2988, DOI: 10.1093/bioinformatics/bts545 (2012).
  * [9] Peters, T. J. _et al._ De novo identification of differentially methylated regions in the human genome.  _Epigenetics & chromatin_ 8, 6, DOI: 10.1186/1756-8935-8-6 (2015).
  * [10] Jaffe, A. E. _et al._ Bump hunting to identify differentially methylated regions in epigenetic epidemiology studies.  _International journal of epidemiology_ 41, 200–209, DOI: 10.1093/ije/dyr238 (2012).
  * [11] Krause, L. _et al._ Identification of the cimp-like subtype and aberrant methylation of members of the chromosomal segregation and spindle assembly pathways in esophageal adenocarcinoma.  _Carcinogenesis_ 37, 356–365, DOI: 10.1093/carcin/bgw018 (2016).
  * [12] Yu, M. _et al._ Subtypes of barrett’s oesophagus and oesophageal adenocarcinoma based on genome-wide methylation analysis.  _Gut_ 68, 389–399, DOI: 10.1136/gutjnl-2017-314544 (2019).
  * [13] Brown, M. B. & Forsythe, A. B.  Robust tests for equality of variances.  _J Am Stat Assoc_ 69, 364–367 (1974).
  * [14] Moinova, H. R. _et al._ Identifying dna methylation biomarkers for non-endoscopic detection of barrett’s esophagus.  _Science translational medicine_ 10, DOI: 10.1126/scitranslmed.aao5848 (2018).
  * [15] Nan, P. _et al._ Mta1 promotes tumorigenesis and development of esophageal squamous cell carcinoma via activating the mek/erk/p90rsk signaling pathway.  _Carcinogenesis_ 41, 1263–1272, DOI: 10.1093/carcin/bgz200 (2020).
  * [16] Toh, Y. _et al._ Overexpression of the mta1 gene in gastrointestinal carcinomas: correlation with invasion and metastasis.  _International journal of cancer_ 74, 459–463, DOI: 10.1002/(sici)1097-0215(19970822)74:4<459::aid-ijc18>3.0.co;2-4 (1997).
  * [17] Zhang, Y. _et al._ Circdido1 inhibits gastric cancer progression by encoding a novel dido1-529aa protein and regulating prdx2 protein stability.  _Molecular cancer_ 20, 101, DOI: 10.1186/s12943-021-01390-y (2021).
  * [18] Peters, T. J. _et al._ Calling differentially methylated regions from whole genome bisulphite sequencing with dmrcate.  _Nucleic acids research_ 49, e109, DOI: 10.1093/nar/gkab637 (2021).