# Adaptive Bayesian Variable Clustering via Structural Learning of Breast
Cancer Data

Elsevier111Since 1880. Radarweg 29, Amsterdam  Riddhi Pratim Ghosh
<EMAIL_ADDRESS>Arnab Kumar Maity
<EMAIL_ADDRESS>Mohsen Pourahmadi<EMAIL_ADDRESS>Bani K.
Mallick<EMAIL_ADDRESS>Department of Biostatistics, University of
Pennsylvania  Early Clinical Development Oncology Statistics, Pfizer, 10777
Science Center Drive, San Diego, California 92121  Department of Statistics,
Texas A&M University, 3143 TAMU, College Station, Texas 77843

###### Abstract

Clustering of proteins is of interest in cancer cell biology. This article
proposes a hierarchical Bayesian model for protein (variable) clustering
hinging on correlation structure. Starting from a multivariate normal
likelihood, we enforce the clustering through prior modeling using angle based
unconstrained reparameterization of correlations and assume a truncated
Poisson distribution (to penalize the large number of clusters) as prior on
the number of clusters. The posterior distributions of the parameters are not
in explicit form and we use a reversible jump Markov chain Monte Carlo
(RJMCMC) based technique is used to simulate the parameters from the
posteriors. The end products of the proposed method are estimated cluster
configuration of the proteins (variables) along with the number of clusters.
The Bayesian method is flexible enough to cluster the proteins as well as the
estimate the number of clusters. The performance of the proposed method has
been substantiated with extensive simulation studies and one protein
expression data with a hereditary disposition in breast cancer where the
proteins are coming from different pathways.

###### keywords:

Angular Reparameterization, Bayesian clustering, Pathways, Reversible jump
Markov chain Monte Carlo

††journal: Journal of LaTeX Templatesmytitlenotemytitlenotefootnotetext: Fully
documented templates are available in the elsarticle package on CTAN.

##  1 Introduction

In cell biology different pathways emerge as they play different and critical
role in cell functions. Even though the functionality of a cell is an outcome
of all the pathway protein expressions as a whole, the individual analysis of
each protein has the potential to unveil the complex characterization of the
cell biology which is the key to understand the proper cell function [Ben-Dor
et al., 1999] .

The goal of the clustering is to distill the data down to a more
comprehensible level subdividing the omics data [D’haeseleer, 2005]. In this
article we focus on the proteomics data and the interest is to infer about the
pathways based on the proteins data via the clustering technique. Clustering
of proteins is a form of unsupervised learning where the proteins are grouped
on the basis of some similarity measures inherent among them. Such clusters
can be mapped to find the appropriate pathway based on the available protein
expression data. Given the functions of the proteins which are measured via
the RPPA technology based protein expressions, it is of interest to track back
the pathways in which the group of proteins belong to, assuming that the
pathways do not have an overlap.

A proper clustering method which explicates the pattern involved in the gene
expression depending on the over-expression or under-expression of those
uncovers the tumor subtypes. For example, Pollack et al. [2002] showed that a
profiling of DNA copy number variation has the potential to detect more
aggressive breast tumors. Washburn et al. [2003] considered correlation of
mRNA and protein expression of amino acid and nucleotide biosynthetic pathway
components for clustering. Ben-Dor et al. [1999] provided an algorithm Cluster
Affinity Search Technique which uses an affinity measure between nodes of a
graph where the genes are represented as the nodes of a graph. In the absence
of genuine variable clustering methods, very often traditional data clustering
algorithms have been applied to this setup using brute force [Vigneau and
Qannari, 2003, Duda et al., 2001] or ad-hoc algorithms based on aspects of
correlation matrices have been proposed. We refer the readers to Jiang et al.
[2004] for a detailed discussion of various correlation based clustering
approaches which have been previously used in literature for analyzing
different gene expression data. However, in the current era of next generation
sequencing the amount of data that one receives and underlying complexity of
the pattern often pose challenges for interpretation and understanding the
results, necessitating a proper and meaningful clustering tool.

In this article, our aim is to cluster the proteins, essentially a variable
clustering technique which is drastically different from approaches for
clustering observations or subjects. To understand it better, let ${\bf Y}$
denote a $n\times k$ data matrix consisting of $k$ proteins and $n$ patients,
represented in the matrix form

$\bf Y=\bordermatrix{&\text{Protein 1}&\text{Protein 2}&\text{Protein
3}&\text{Protein 4}&\dots&\text{Protein
k}\cr&y_{11}&y_{12}&y_{13}&y_{14}&\dots&y_{1k}\cr&y_{21}&y_{22}&y_{23}&y_{24}&\dots&y_{2k}\cr&y_{31}&y_{32}&y_{33}&y_{34}&\dots&y_{3k}\cr&\vdots&\vdots&\vdots&\ddots&\vdots\cr&y_{n1}&y_{n2}&y_{n3}&y_{n4}&\dots&y_{nk}}$
(1)

From $\eqref{data_matrix_eqn}$, one notes that each of $n$ rows corresponds to
one patient and each of $k$ columns pertains to one protein. A typical data
clustering approach partitions the rows of $\bf Y$, i.e. essentially
clustering of the patients. We are interested in partitioning the columns of
$\bf Y$ which is essentially clustering of the proteins, and correlations
between the proteins serve as our main building block to implement the
algorithm. In a typical data clustering algorithm we consider how similar the
objects are based on a similarity norm (say Euclidean or some other kind of
distance). On the contrary, in a variable clustering problem, we are concerned
with the correlation among the variables. Hence, highly correlated proteins
are more likely to lie in the same cluster. As an example, consider a cluster
analysis of a set of proteins which belong to different signaling pathways
assuming the pathways are not overlapping. The genetic behaviors control the
proliferation of a cell or death of a cell; and depending on signals the
proteins receive and send, the cell structure is classified into signaling
pathways. In turn, one can assume that the similarly expressed proteins belong
to the same pathway which can be recovered via a variable cluster analysis.

Among the different algorithmic clustering techniques commonly used in
practice, hierarchical clustering (agglomerative and divisive approach) and
partition methods (K-means clustering) hing on a distance metric [Bibby et
al., 1979, Friedman et al., 2001, Rokach and Maimon, 2005] without assuming
any underlying probability model for the clusters. In addition, model based
approach usually assumes a mixture model for the data. Even though there is a
vast amount of works in the field of data clustering, but the variable
clustering problem is at its infancy and has gotten limited attention [Bunea
et al., 2020]. The literature on Bayesian methods for variable clustering is
also sparse with a few notable exceptions (Liechty et al. 2004, Palla et al.
2012). Palla et al. [2012] developed a nonparametric Bayes algorithm based on
Chinese restaurant process. On the other hand, our method is in the spirit of
Liechty et al. [2004] where a parametric model based approach has been
considered. A key advantage of our approach is that the number of clusters is
assumed to be unknown apriori, and is determined using a reversible jump
Markov Chain Monte Carlo algorithm [RJMCMC] [Green, 1995].

In this article, our contributions can be summarized as, first, to develop the
model-based variable clustering method with block common correlation
structures. Second, we propose a novel variable clustering algorithm using the
angular representation of the correlations [Pinheiro and Bates, 1996,
Rapisarda et al., 2007, Tsay and Pourahmadi, 2017, Ghosh et al., 2020] and the
ensuing angles (hyperspherical coordinates). Third, we elicit substantive
prior information on these angles which makes clustering of the variables
feasible, a data-driven estimate of number of clusters which traditional
algorithms fail to provide. For the posterior inference, since the angle
parameters are badly entangled in the posterior distribution, we resort to the
Markov chain Monte Carlo algorithm [Tierney, 1994]. For the posterior
inference, we resort to the standard RJMCMC techniques as in [Green, 1995,
Robert, 2004, Green and Hastie, 2009, Fan and Sisson, 2011]. The rest of the
article is organized as follows. In 2 we review angular reparameterization of
a correlation matrix and present clustering model through prior specification
on the angles. In section 3 we describe our posterior computation through
RJMCMC. Section 4 presents simulation results and clustering of a protein
expression data. Finally section 5 concludes the article.

##  2 Review of angular reparametrization ($\Theta$) of $R$

This section describes connections between the hyperspherical coordinates
(angles) and a correlation matrix $R=(r_{ij})$.

For a general $k\times k$ correlation matrix $R$ with $1$’s in the diagonal,
its Cholesky decomposition is given by $R=BB^{\top}$ where the Cholesky factor
$B$ is a lower triangular matrix. Since the rows of $B$ are vectors of unit-
length, it turns out that they admit the following representation involving
trigonometric functions of some angles [Pinheiro and Bates, 1996, Rapisarda et
al., 2007]:

$B=\begin{bmatrix}1&0&0&0&\dots&0\\\ c_{21}&s_{21}&0&0&\dots&0\\\
c_{31}&c_{32}s_{31}&s_{32}s_{31}&0&\dots&0\\\
c_{41}&c_{42}s_{41}&c_{43}s_{42}s_{41}&\prod_{j=1}^{3}s_{4j}&\dots&0\\\
\vdots&\vdots&\vdots&\ddots&\vdots\\\
c_{k1}&c_{k2}s_{k1}&c_{k3}s_{k2}s_{k1}&c_{k4}\prod_{j=1}^{3}s_{kj}&\dots&\prod_{j=1}^{k-1}s_{kj}\\\
\end{bmatrix}$ (2)

with $c_{ij}=\text{cos}(\theta_{ij})$ and $s_{ij}=\text{sin}(\theta_{ij})$,
where the angles $\theta_{ij}$’s are measured in radians, $1\leq j<i\leq k$.
Restricting $\theta_{ij}\in[0,\pi)$ makes the diagonal entries of $B$ non-
negative, and hence $B$ is unique to which we associate a $(k-1)\times(k-1)$
lower triangular matrix $\Theta$ with $k(k-1)/2$ angles:

$\Theta=\begin{bmatrix}\theta_{21}&0&0&\dots&0\\\
\theta_{31}&\theta_{32}&0&\dots&0\\\ \vdots&\vdots&\vdots&\ddots&\vdots\\\
\theta_{k1}&\theta_{k2}&\theta_{k3}&\dots&\theta_{k,k-1}\end{bmatrix}$ (3)

Note that the $(i,j)$-th element of $\Theta$ is denoted by $\theta_{i+1,j}$ so
that $\theta_{ij}$ corresponds to the $(i,j)$-th element of $R$, we refer to
$\Theta$ as the angular matrix associated to $R$. For further details and
applications of these angles, see Creal et al. [2011], Zhang et al. [2015],
Tsay and Pourahmadi [2017], and Ghosh et al. [2020] . One can characterize
block diagonal correlation matrices in terms of structured $\Theta$ matrix,
which is completely determined by some (pivotal) angles.

###  2.1 Correspondence of clustering between $R$ and $\Theta$

######  Proposition 2.1.

For a block diagonal correlation matrix $R=\text{block
diag}(R_{1},R_{2},\cdots,R_{m})$, consisting of $m$ equicorrelated blocks
($r_{i}$ for block $R_{i}$), the corresponding angular matrix $\Theta$ is
characterized by only $m$ angles $\theta_{1},\theta_{2},...,\theta_{m}$, where
$r_{i}=\text{cos }\theta_{i}$.

###### Proof.

See section A.1. ∎

It follows immediately from proposition 2.1 that in case of block diagonal
correlation matrix, clustering on correlations rendering to $m$ different
groups is equivalent to clustering of those $m$ angles by the monotonicity of
cosine function. However, this will be impose some conditions on the pivotal
angles to maintain positive definiteness. Assuming each block has dimension
$k_{i}$ so that $\sum_{i=1}^{m}k_{i}=k$, the support of $\theta_{i}$ is
$0<\theta_{i}<\text{arccos}(1/(k_{i}-1))$ for $i=1,2,\cdots,m$.

######  Proposition 2.2.

Suppose that $r_{1}=\text{cos }\theta_{1}$, $r_{2}=\text{cos }\theta_{2}$.
Then $|\theta_{1}-\theta_{2}|\geq\delta$ if and only if
$|r_{1}-r_{2}|\geq|1-\text{cos }\delta|$

###### Proof.

See section A.2. ∎

###  2.2 Likelihood function

In this article, we assume throughout that the data ${\bf y_{1}},{\bf
y_{2}},\cdots,{\bf y_{n}}$ follow a zero mean normal distribution with
covariance assumed to be the correlation matrix $R$. As noted in proposition
2.1, $R$ can be written as a function of $m$ pivotal angles,
$\theta_{pivot}=(\theta_{1},\theta_{2},...,\theta_{m})^{\top}$, and hence
denoting the transformation from $\theta_{pivotal}$ to $R$ by $T$, the
likelihood is proportional to,

$\displaystyle L({\bf y_{1}},{\bf y_{2}},...,{\bf
y_{n}}|\theta_{pivot})\propto\text{ det }(T(\theta_{pivot}))^{-n/2}\text{ exp
}\\{-\frac{1}{2}ST^{-1}(\theta_{pivot})\\}$ (4)

where $S=\sum_{i=1}^{n}{\bf y_{i}}{\bf y_{i}}^{\top}$.

###  2.3 Prior specification on the angles

The number of clusters m can take any value any value in 1,2,…,k. Therefore,
we assume a truncated Poisson distribution on m. Given m, define $k\times m$
matrix Z whose i-th row corresponds to the allocation of $i$-th variable in
one of the $m$ clusters, i.e.

${\bf Z}_{iu}=\begin{cases}1\;\;\;\ \text{if $i$-th variable belongs to $u$-th
cluster}\\\ 0\;\;\ \text{otherwise }\end{cases}$

Since we are assuming that a variable belongs to exactly one cluster,
therefore, each row of ${\bf Z}$ contains exactly one 1 and rests are 0s. We
assume the following hierarchical prior models for model parameters by
assuming a truncated Poisson distribution on number of clusters $m$ to
penalize large number of clusters, multinomial distribution on each row of the
indicator matrix ${\bf Z}$ and a Dirichlet distribution for the multinomial
hyper-parameters. The hierarchical prior structure is succinctly described as,

$\displaystyle m\sim truncPois(m;1,k)$ (5) $\displaystyle{\bf
q}=(q_{1},q_{2},...,q_{m})^{\top}\sim
Dirichlet(\alpha_{1},\alpha_{2},...,\alpha_{m})$ (6) $\displaystyle{\bf
Z}_{i}\sim Multinomial(1;q_{1},q_{2},...,q_{m})\;\;\ \text{ for }i=1,2,...,k,$
(7)

where $\alpha_{i}$s are any positive numbers and $truncPois(m;1,k)$ is a
truncated Poisson distribution supported on the integers in between 1 and $k$
(number of proteins or variables) for the number of clusters $m$. Having
sampled ${\bf Z}$, the allocations are determined. Let $k_{u}$ denote the size
of $u$-th cluster,

$\displaystyle k_{u}=|\\{i:z_{iu}=1\\}|$ (8)

for $u=1,2,\cdots,m$.

Then assume the following prior on
$\theta_{piv}=(\theta_{1},\theta_{2},...,\theta_{m})^{\top}$ in order to
shrink them to different values.

$\displaystyle\theta_{piv}|{\bf
Z},m,\Lambda=\prod_{u=1}^{m}Q\Big{(}\theta_{u};0,\text{arccos}\big{(}\frac{1}{k_{u}-1}\big{)},\lambda_{u}\Big{)}$
(9)

where $Q(\theta;0,a,\lambda)$ is the density of truncated wrapped Exponential
distribution [Mardia and Jupp, 2009] between 0 and $a$ with parameter
$\lambda$. We are clustering the pivotal angles by introducing wrapped
exponential distribution distribution with different parameters. Suppose
$\Lambda=(\lambda_{1},\lambda_{2},...,\lambda_{m})^{\top}$ and we sample
$\lambda_{1},\lambda_{2},\cdots,\lambda_{m}$ in the following manner,

$\displaystyle\lambda_{1}\sim N^{+}(\lambda;0,1,0,\infty)$ (10)
$\displaystyle\lambda_{2}|\lambda_{1}\sim
N^{+}(\lambda;0,1,\lambda_{1},\infty)$
$\displaystyle\lambda_{i}|\lambda_{i-1}\sim
N^{+}(\lambda;0,1,\lambda_{i-1},\infty)\;\;\;\text{ for }i=2,3,\cdots,m$ ,

where $N^{+}(;0,1,a,\infty)$ denotes a truncated normal distribution on
$(0,\infty)$ with mean 0 and variance 1 which has the following density,

$\displaystyle
f(\lambda;\mu=0,\sigma=1,a,\infty)=\frac{\phi(\lambda)}{1-\Phi(a)}\;\;\;\ ,$

where $\phi$ and $\Phi$ are the density and distribution function of a
standard normal distribution respectively.

The salient features of the prior formulation of $\lambda_{i}$s given in (10)
are the followings: (1) The prior mean for the $i$-th pivotal angles is
$\mathbb{E}\theta_{i}=\text{arctan}(1/\lambda_{i})\;\;\;\ \text{for
}i=1,2,...,m$. Since these angles vary in $[0,\pi)$, $\lambda_{i}$s take value
on positive real line. (2) Also $\lambda_{i}$’s satisfy
$\lambda_{1}<\lambda_{2}<\lambda_{2}<...<\lambda_{m}$, which enforces
separation of clusters through prior model.

##  3 Posterior computation

With the likelihoood function (4) and prior specified in 2.3, the posterior
distribution is proportional to

$\displaystyle p(\Theta,{\bf Z},\Lambda,m|{\bf y_{1}},{\bf y_{2}},...,{\bf
y_{n}})$ (11) $\displaystyle\propto L({\bf y_{1}},{\bf y_{2}},...,{\bf
y_{n}}|\Theta,{\bf Z},m)\times p(m)\times p({\bf q}|m)\times p({\bf Z}|{\bf
q},m)$ $\displaystyle\times p(\Theta|\Lambda,{\bf Z},m)\times p(\Lambda)$

Our goal in this section is to estimate number of clusters $m$ and posterior
of ${\bf Z}$. The algorithm is, thus, accomplished by performing a reversible
jump Markov chain Monte Carlo(RJMCMC) algorithm.

From proposed priors, one can note that the clusters are induced by the
elements of $\Lambda$, thus, in the following RJMCMC algorithm [Green, 1995,
Robert, 2004, Green and Hastie, 2009, Fan and Sisson, 2011], at each iteration
either one element of $\Lambda$, say $\lambda_{j}$ is randomly split into
$(\lambda_{j_{1}},\lambda_{j_{2}})$ (Birth step) or two elements of $\Lambda$
are merged into a single element (Death step). The algorithm is summarized as
follows.

  * Step 1. Initialize $\Theta$, $\Lambda$. In the initialization step, one may assume any block diagonal correlation structure to initialize $\Theta$.

  * Step 2. A particular iteration, say $q$-th iteration consists of a Birth step and a Death step.

    * Birth Step: Split $\lambda^{(q)}_{j}$ to $(\lambda^{(q)}_{j_{1}},\lambda^{(q)}_{j_{2}})^{\top}$ by $\lambda^{(q)}_{j_{1}}=\lambda^{(q)}_{j}+\tau$, $\lambda^{(q)}_{j_{2}}=\lambda^{(q)}_{j}-\tau$, where $\tau\sim Unif(-\pi/4,\pi/4)$ and dimension of $\lambda^{(q)}$ is increased by 1 with acceptance probability $\alpha=\text{min}\\{1,\frac{p(\Theta^{(q)},\lambda^{(q)}_{j_{1}},\lambda^{(q)}_{j_{2}},d(j_{1},j_{2}))}{p(\Theta^{(q)},\lambda^{(q)}_{j},d(j))}\times\frac{2}{\pi}\times|\frac{\partial(\lambda^{(q)}_{j_{1}},\lambda^{(q)}_{j_{2}})}{\partial(\lambda^{(q)}_{j},\tau)}|\\}$

    * Death step: Two components $\lambda^{(q)}_{j_{1}}$ and $\lambda^{(q)}_{j_{2}}$ are merged to a single component $\lambda^{(q)}_{j}=(\lambda^{(q)}_{j_{1}}-\tau+\lambda^{(q)}_{j_{2}}+\tau)/2$ with acceptance probability $\alpha=\text{min}\\{1,\frac{p(\Theta^{(q)},\lambda^{q}_{j},d(j))}{p(\Theta^{(q)},\lambda^{(q)}_{j_{1}},\lambda^{(q)}_{j_{2}},d(j_{1},j_{2}))}\times\frac{\pi}{2}\times|\frac{\partial(\lambda^{(q)}_{j},\tau)}{\partial(\lambda^{(q)}_{j_{1}},\lambda^{(q)}_{j_{2}})}|\\}$

  * Step 4. Step 1, 2 and 3 are repeated as many times as required to ensure convergence and the value of $m$ is determined by which stage is visited maximum number of times, maximum aposteriori estimate (MAP) and posterior estimate of ${\bf Z}$ is obtained by averaging over those stages.

##  4 Simulations and Data Analyses

In this section, we compare numerical performance of our Bayesian Variable
Clustering (BVC) algorithm with a recent method based on COD (Covariance
Difference) of Bunea et al. [2020], Partitioning Around Medoids (PAM)
algorithm which minimizes the Manhattan distance of the data points to the
medoids [Kaufman and Rousseeuw, 2009] and the classical or standard K-means
clustering algorithm. The performance criterion we use is the proportion of
true recovery which is defined for a $k$ variable as

$\frac{\\#\text{ of variables in the true clusters}}{k}.$ (12)

COD and PAM have been implemented using the R packages cord [Luo et al., 2015]
and class [Venables and Ripley, 2002] available via CRAN and K-means algorithm
has been implemented on the transposed data matrix using kmeans() function in
R software [R Core Team, 2020]. It is instructive to note that quantity in
(12) takes value in the interval [0,1]. As the value approaches to 1, the
recovery becomes better.

###  4.1 Simulation study

We start with an $m\times m$ matrix $C=B^{\top}B$ where the entries of the
random $(m-1)\times m$ matrix $B$ take values $-1,0,1$ with probabilities
$0.5\times m^{-1/2}$, $1-m^{-1/2}$ and $0.5\times m^{-1/2}$, respectively,
with $m$ being the number of clusters. Next, we consider a balanced case with
each group (cluster) of size $k/m$. Let $A=(a_{ij})$ be the $k\times m$
membership matrix with $a_{ij}=1$ if the $i$-th variable belongs to $C_{j}$
and 0 otherwise. Finally, consider the covariance matrix
$\Sigma=ACA^{\top}+\Gamma$ where $\Gamma$ is a diagonal matrix whose entries
are random permutations of $\\{0.5,0.5+1.5/(k-1),...,2\\}$ and the
corresponding correlation matrix $R$. With $k=200,m=4$, we simulate $n$
independent observations from a multivariate normal distribution with mean
zero vector and covariance matrix $R$, where we vary $n$ in 100, 300, 600, 900
to compare BVC, COD and K-means algorithms with respect to cluster recovery
criterion in (12). The results presented in Figure 1 shows the superior
performance of BVC relative to COD, PAM and K-means as the values stay closer
to 1.

Figure 1: Comparing BVC(blue), COD(red), PAM(green) and K-means(black) for
simulation study in section 4.1

###  4.2 Application of Protein clustering to Hereditary Breast Cancer Data

Breast cancer is one of the most common cancers with a massive number of cases
reported. For instance, in 2018, more than 268,000 Americans were estimated to
have been diagnosed and 41,000 were estimated to have died from breast cancer
related tumors [Bray et al., ]. The Cancer Genome Atlas: TCGA is the largest
available cancer data consortium consisting of parallel mRNA expressions, DNA
copy number, methylation expressions, protein expressions, along with clinical
variables such as survival or the tumor stages for a total of 33 types of
tumors. Among them we consider the information of 222 breast tumor samples; we
consider 27 different proteins 4 different pathways (see section A.3).
Different subsets of this data has been used in Maity et al. [2020a, b] for
different purposes.

Applying our BVC algorithm to this data, the MAP estimate of the number of
clusters is 4, which is consistent with the number of pathways. However,
applying the COD algorithm in Bunea et al. [2020] the estimated number of
clusters is 23, much larger than the known value of 4. In section A.4, we
provide the assignments of various proteins in different clusters.
Additionally, for the sake of comparison we have also applied the K-means
algorithm to this data for $k=4,23$, respectively, with results reported in
the section A.4. The results suggest that our Bayesian variable clustering
(BVC) is performing better to cluster the proteins with respect to pathways.
Only misclassified proteins are MAPK_pT201_Y204, CD31, CD49b, CDK1. The COD
algorithm reports that number of clusters is 23 which appears to be too high
since the number of proteins is 27. A possible reason could be this algorithm
is meant for high dimensional clustering, it fails to detect clustering
configuration in small dimensional cases. Comparisons with standard K-means
and PAM algorithm also reveal that these two methods result in more
disagreement of the cluster configuration of the proteins according to the
pathway information. This apart, K-means and PAM algorithm disagree among
themselves, e.g., ER-alpha, JNK_pT183_pT185 etc.(Table 2). We have also
performed hierarchical clustering on this data with various linkages . The
results are presented in Figure 2.

Figure 2: Hierarchical clustering for the protein expression data with four
different linkages.

##  5 Discussion

We have proposed a correlation matrix based Bayesian clustering technique to
recover the protein signaling pathways. This method uses angular
reparameterization of correlation matrix with the specification of wrapped
exponential prior on the angle parameters. Nonetheless, as an alternative, one
can use any truncated circular distribution as prior for pivotal angles, for
example von-Mises distribution. However, this particular choice produces a
mean which has no closed form and as a result our proposed method can not be
carried out for a posterior analysis.

A large amount of recent interest is being channelized to analyze the
proteomics data directly because direct analysis of proteins has potential to
uncover the cell functional characteristics. When it is of interest to find
the group of proteins having similar functions which may be evident via their
expression measurements then our proposed method can be used to bridge that
gap. As mentioned earlier, our method is particularly useful when the number
of clusters is not known and hence is learned via the posterior MCMC, which is
often the case for the real data where the determining the number of clusters
is itself a tedious job.

## Acknowledgment

The research reported in this paper was supported by grant from the National
Institutes of Health (R01CA194391).

##  Appendix A Appendix

###  A.1 Proof of proposition 2.1

###### Proof.

The proof uses the Cholesky decomposition of $R$ and the fact that Cholesky
factor of a block diagonal matrix is also block diagonal and vice versa.

Hence, the lower triangular Cholesky factor $R$ has the form

$B=\text{block diag}(B_{1},B_{2},...,B_{m})$, where $B_{i}$ is an upper
triangular Cholesky factor of $R_{i}$.

The proof will be complete if we can show that Cholesky factor $B$ of a
compound symmetric correlation matrix $R$ can be written in terms of only one
angle. From the relationship between angles and Cholesky factor, it follows
that $\text{cos}(\theta_{i1})=b_{i1}=r$ for $j=2,3,...,k$. Thus
$\theta_{i1}$’s are all equal to a common $\theta$. For any $i>j$, the proof
follows by induction. For $i>j$,
$r_{ij}=\sum_{l=1}^{j-1}b_{il}b_{jl}+b_{jj}b_{ij}$. By induction hypothesis,
all the preceding angles and Cholesky factors are functions of $r$. Thus, the
first term is a function of $r$. For the second term, we note
$b_{jj}=\prod_{l=1}^{j-1}\text{sin}(\theta_{jl})$, which involves all the
preceding angles and thus a function of $r$ and
$b_{ij}=\text{cos}(\theta_{ij})\prod_{l=1}^{j-1}\text{sin}(\theta_{jl})$. Thus
it follows that $\text{cos}(\theta_{ij})$ is a function of $r$ and
$\theta_{ij}$ is a function of $\theta$.

∎

###  A.2 Proof of proposition 2.2

###### Proof.

First consider $|\theta_{1}-\theta_{2}|=\delta$. Note that
$|r_{1}-r_{2}|=|\int_{\theta_{1}}^{\theta_{2}}\text{sin }xdx|$. Also it is
clear that $|r_{1}-r_{2}|$ is an increasing function of
$|\theta_{1}-\theta_{2}|$, since sin is positive in $[0,\pi)$. Now since sin
is increasing in $[0,\pi/2]$ and decreasing in $(\pi/2,\pi)$, $|r_{1}-r_{2}|$
will take minimum value for $|\theta_{1}-\theta_{2}|=\delta$ when
$\theta_{1}=0,\theta_{2}=\delta$. Thus the minimum value of $|r_{1}-r_{2}|$ is
$|1-\text{cos}\delta|$. ∎

###  A.3 Pathway Information

See Table 1 for pathway information of the proteins.

Table 1: Pathway Protein List MAP kinase | PI3K/AKT/mTOR | JAK-STAT | Wnt
---|---|---|---
ER-alpha | AKT | SHC_pY317 | CD31
ER-alpha_pS118 | AKT_pS473 | STAT3_pY705 | CD49b
ERK2 | AKT_pT308 | STAT5-alpha | CDK1
JNK2 | FOXO3a |  | Cyclin_D1
JNK_pT183_pT185 | PTEN |  | Fibronectin
MAPK_pT202_Y204 | mTOR |  | GSK3-alpha-beta
p38_MAPK | mTOR_pS2448 |  | GSK3-alpha-beta_pS21_S9
p38_pT180_Y182 |  |  | VEGFR2
|  |  | beta-Catenin

###  A.4 Cluster Assignments of proteins

Table 2 presents the cluster assignments of proteins by BVC, COD, PAM and
K-means algorithms.

Table 2: Cluster comparisons by BVC, COD and K-means Protein | BVC | COD | K-means ($m$=4) | PAM ($m$=4) | K-means (k=23)
---|---|---|---|---|---
ER-alpha | $C_{1}$ | $C_{1}$ | $C_{3}$ | $C_{1}$ | $C_{2}$
ER-alpha_pS118 | $C_{1}$ | $C_{1}$ | $C_{1}$ | $C_{2}$ | $C_{7}$
ERK2 | $C_{1}$ | $C_{2}$ | $C_{1}$ | $C_{2}$ | $C_{15}$
JNK2 | $C_{1}$ | $C_{5}$ | $C_{1}$ | $C_{2}$ | $C_{10}$
JNK_pT183_pT185 | $C_{1}$ | $C_{6}$ | $C_{4}$ | $C_{2}$ | $C_{8}$
MAPK_pT202_Y204 | $C_{3}$ | $C_{7}$ | $C_{2}$ | $C_{3}$ | $C_{112}$
p38_MAPK | $C_{1}$ | $C_{8}$ | $C_{4}$ | $C_{2}$ | $C_{9}$
p38_pT180_Y182 | $C_{1}$ | $C_{9}$ | $C_{4}$ | $C_{2}$ | $C_{22}$
AKT | $C_{2}$ | $C_{10}$ | $C_{1}$ | $C_{2}$ | $C_{11}$
AKT_pS473 | $C_{2}$ | $C_{11}$ | $C_{2}$ | $C_{4}$ | $C_{1}$
AKT_pT308 | $C_{2}$ | $C_{12}$ | $C_{2}$ | $C_{4}$ | $C_{1}$
FOXO3a | $C_{2}$ | $C_{13}$ | $C_{4}$ | $C_{2}$ | $C_{5}$
PTEN | $C_{2}$ | $C_{14}$ | $C_{1}$ | $C_{2}$ | $C_{13}$
mTOR | $C_{2}$ | $C_{15}$ | $C_{1}$ | $C_{2}$ | $C_{21}$
mTOR_pS2448 | $C_{2}$ | $C_{16}$ | $C_{1}$ | $C_{2}$ | $C_{20}$
SHC_pY317 | $C_{3}$ | $C_{17}$ | $C_{4}$ | $C_{2}$ | $C_{17}$
STAT3_pY705 | $C_{3}$ | $C_{18}$ | $C_{4}$ | $C_{2}$ | $C_{23}$
STAT5-alpha | $C_{3}$ | $C_{2}$ | $C_{1}$ | $C_{2}$ | $C_{18}$
CD31 | $C_{2}$ | $C_{3}$ | $C_{4}$ | $C_{2}$ | $C_{6}$
CD49b | $C_{2}$ | $C_{19}$ | $C_{4}$ | $C_{2}$ | $C_{5}$
CDK1 | $C_{2}$ | $C_{3}$ | $C_{4}$ | $C_{2}$ | $C_{5}$
Cyclin_D1 | $C_{4}$ | $C_{20}$ | $C_{4}$ | $C_{2}$ | $C_{16}$
Fibronectin | $C_{4}$ | $C_{21}$ | $C_{4}$ | $C_{2}$ | $C_{4}$
GSK3-alpha-beta | $C_{4}$ | $C_{4}$ | $C_{1}$ | $C_{2}$ | $C_{21}$
GSK3-alpha-beta_pS21_S9 | $C_{4}$ | $C_{22}$ | $C_{2}$ | $C_{4}$ | $C_{19}$
VEGFR2 | $C_{4}$ | $C_{23}$ | $C_{1}$ | $C_{2}$ | $C_{3}$
beta-Catenin | $C_{4}$ | $C_{4}$ | $C_{1}$ | $C_{2}$ | $C_{14}$

## References

  * Ben-Dor et al. [1999] Ben-Dor, A., Shamir, R., Yakhini, Z., 1999.  Clustering gene expression patterns.  Journal of Computational Biology 6, 281–297.
  * Bibby et al. [1979] Bibby, J., Kent, J., Mardia, K., 1979.  Multivariate Analysis.
  * [3] Bray, F., Ferlay, J., Soerjomataram, I., Siegel, R.L., Torre, L.A., Jemal, A., .  Global Cancer statistics 2018: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries .
  * Bunea et al. [2020] Bunea, F., Giraud, C., Luo, X., Royer, M., Verzelen, N., et al., 2020.  Model assisted variable clustering: minimax-optimal recovery and algorithms.  The Annals of Statistics 48, 111–137.
  * Creal et al. [2011] Creal, D., Koopman, S.J., Lucas, A., 2011.  A dynamic multivariate heavy-tailed model for time-varying volatilities and correlations.  Journal of Business & Economic Statistics 29, 552–563.
  * D’haeseleer [2005] D’haeseleer, P., 2005.  How does gene expression clustering work?  Nature Biotechnology 23, 1499\.
  * Duda et al. [2001] Duda, R.O., Hart, P.E., Stork, D.G., 2001.  Pattern Classification.
  * Fan and Sisson [2011] Fan, Y., Sisson, S.A., 2011\.  Reversible jump MCMC.  Handbook of Markov Chain Monte Carlo , 67–92.
  * Friedman et al. [2001] Friedman, J., Hastie, T., Tibshirani, R., 2001.  The Elements of Statistical Learning. volume 1.  Springer Series in Statistics New York.
  * Ghosh et al. [2020] Ghosh, R.P., Mallick, B.K., Pourahmadi, M., 2020.  Bayesian Estimation of Correlation Matrices of Longitudinal Data.  Bayesian Analysis .
  * Green [1995] Green, P.J., 1995.  Reversible jump Markov chain Monte Carlo computation and Bayesian model determination.  Biometrika 82, 711–732.
  * Green and Hastie [2009] Green, P.J., Hastie, D.I., 2009\.  Reversible jump MCMC.  Genetics 155, 1391–1403.
  * Jiang et al. [2004] Jiang, D., Tang, C., Zhang, A., 2004.  Cluster Analysis for Gene Expression Data: A Survey.  IEEE Transactions on Knowledge and Data Engineering 16, 1370–1386.
  * Kaufman and Rousseeuw [2009] Kaufman, L., Rousseeuw, P.J., 2009\.  Finding Groups in Data: an Introduction to Cluster Analysis. volume 344.  John Wiley & Sons.
  * Liechty et al. [2004] Liechty, J.C., Liechty, M.W., Müller, P., 2004.  Bayesian correlation estimation.  Biometrika 91, 1–14.
  * Luo et al. [2015] Luo, X., Bunea, F., Giraud, C., 2015.  cord: Community Estimation in G-Models via CORD.  URL: https://CRAN.R-project.org/package=cord. r package version 0.1.1.
  * Maity et al. [2020a] Maity, A.K., Bhattacharya, A., Mallick, B.K., Baladandayuthapani, V., 2020a.  Bayesian data integration and variable selection for pan-cancer survival prediction using protein expression data.  Biometrics 76, 316–325.
  * Maity et al. [2020b] Maity, A.K., Lee, S.C., Mallick, B.K., Sarkar, T.R., 2020b.  Bayesian Structural Equation Modeling in Multiple Omics Data with Application to Circadian Genes.  Bioinformatics 36, 3951 – 3958.
  * Mardia and Jupp [2009] Mardia, K.V., Jupp, P.E., 2009\.  Directional Statistics. volume 494\.  John Wiley & Sons.
  * Palla et al. [2012] Palla, K., Ghahramani, Z., Knowles, D.A., 2012.  A nonparametric variable clustering model, in: Advances in Neural Information Processing Systems, pp. 2987–2995.
  * Pinheiro and Bates [1996] Pinheiro, J.C., Bates, D.M., 1996\.  Unconstrained parametrizations for variance-covariance matrices.  Statistics and Computing 6, 289–296.
  * Pollack et al. [2002] Pollack, J.R., Sørlie, T., Perou, C.M., Rees, C.A., Jeffrey, S.S., Lonning, P.E., Tibshirani, R., Botstein, D., Børresen-Dale, A.L., Brown, P.O., 2002\.  Microarray analysis reveals a major direct role of DNA copy number alteration in the transcriptional program of human breast tumors.  Proceedings of the National Academy of Sciences 99, 12963–12968.
  * R Core Team [2020] R Core Team, 2020.  R: A Language and Environment for Statistical Computing.  R Foundation for Statistical Computing. Vienna, Austria.  URL: https://www.R-project.org/.
  * Rapisarda et al. [2007] Rapisarda, F., Brigo, D., Mercurio, F., 2007.  Parameterizing correlations: a geometric interpretation.  IMA Journal of Management Mathematics 18, 55–73.
  * Robert [2004] Robert, C.P., 2004.  Monte Carlo methods.  Wiley Online Library.
  * Rokach and Maimon [2005] Rokach, L., Maimon, O., 2005\.  Clustering methods, in: Data Mining and Knowledge Discovery Handbook. Springer, pp. 321–352.
  * Tierney [1994] Tierney, L., 1994.  Markov chains for exploring posterior distributions.  the Annals of Statistics , 1701–1728.
  * Tsay and Pourahmadi [2017] Tsay, R.S., Pourahmadi, M., 2017\.  Modelling structured correlation matrices.  Biometrika 104, 237–242.
  * Venables and Ripley [2002] Venables, W.N., Ripley, B.D., 2002\.  Modern Applied Statistics with S.  Fourth ed., Springer, New York.  URL: http://www.stats.ox.ac.uk/pub/MASS4. iSBN 0-387-95457-0.
  * Vigneau and Qannari [2003] Vigneau, E., Qannari, E., 2003\.  Clustering of variables around latent components.  Communications in Statistics-Simulation and Computation 32, 1131–1150.
  * Washburn et al. [2003] Washburn, M.P., Koller, A., Oshiro, G., Ulaszek, R.R., Plouffe, D., Deciu, C., Winzeler, E., Yates, J.R., 2003\.  Protein pathway and complex clustering of correlated mRNA and protein expression analyses in Saccharomyces cerevisiae.  Proceedings of the National Academy of Sciences 100, 3107–3112.
  * Zhang et al. [2015] Zhang, W., Leng, C., Tang, C.Y., 2015.  A joint modelling approach for longitudinal studies.  Journal of the Royal Statistical Society: Series B (Statistical Methodology) 77, 219–238.