# coseq package: Quick-start guide

#### 2017-04-26

coseq is a package to perform clustering analysis of sequencing data (e.g., co-expression analysis of RNA-seq data), using transformed profiles rather than the raw counts directly. The package implements two distinct strategies in conjunction with transformations: 1) MacQueen’s K-means algorithm and 2) Gaussian mixture models. For both strategies model selection is provided using the slope heuristics approach and integrated completed likelihood (ICL) criterion, respectively. Note that for backwards compatibility, coseq also provides a wrapper for the Poisson mixture model originally proposed in Rau et al. (2015) and implemented in the HTSCluster package.

The methods implemented in this package are described in detail in the following three publications.

• Godichon-Baggioni, A., Maugis-Rabusseau, C. and Rau, A. (2017) Clustering transformed compositional data using K-means, with applications in gene expression and bicycle sharing system data. arXiv:1704.06150. [Paper that introduced the use of the K-means algorithm for RNA-seq profiles after transformation via the centered log ratio (CLR) or log centered log ratio (logCLR) transformation.]

• Rau, A. and Maugis-Rabusseau, C. (2017) Transformation and model choice for RNA-seq co-expression analysis. Briefings in Bioinformatics, doi: http://dx.doi.org/10.1101/065607. [Paper that introduced the idea of clustering profiles for RNA-seq co-expression, and suggested using Gaussian mixture models in conjunction with either the arcsine or logit transformation.]

• Rau, A., Maugis-Rabusseau, C., Martin-Magniette, M.-L., Celeux, G. (2015) Co-expression analysis of high-throughput transcriptome sequencing data with Poisson mixture models. Bioinformatics, 31(9): 1420-1427. link [Paper that introduced the use of Poisson mixture models for RNA-seq counts.]

Below, we provide a quick-start guide using a subset of RNA-seq data to illustrate the functionalities of the coseq package. In this document, we focus on the methods described in Rau and Maugis-Rabusseau (2017) and Godichon-Baggioni et al. (2017). For more information about the method described in Rau et al. (2015), see the HTSCluster vignette.

## Quick start (tl;dr)

A standard coseq analysis takes the following form, where counts represents a matrix or data.frame of gene-level counts arising from an RNA-seq experiment (of dimension n x d for n genes and d samples). Results, exported in the form of a coseqResults S4 object, can easily be examined using standard summary and plot functions (see below and the User’s Guide for examples of plots that may be obtained):

run <- coseq(counts, K=2:25)
summary(run)
plot(run)

The cluster labels obtained from the coseq clustering may easily be obtained using the clusters functions:

clusters(run)

Note that unless otherwise specified by the user, coseq uses the following default parameters:

• K-means algorithm
• Log centered log transformation (logclr) on profiles
• Library size normalization via the Trimmed Mean of M-values (TMM) approach
• Model selection (choice of the number of clusters) via the slope heuristics approach
• No parallelization.

## Co-expression analysis with coseq

For the remainder of this vignette, we make use of the mouse neocortex RNA-seq data from Fietz et al. (2012) (available at https://perso.math.univ-toulouse.fr/maugis/mixstatseq/packages and as an ExpressionSet object called fietz in coseq). The aim in this study was to investigate the expansion of the neocortex in five embryonic (day 14.5) mice by analyzing the transcriptome of three regions: the ventricular zone (VZ), subventricular zone (SVZ) and cortical place (CP).

We begin by first loading the necessary R packages as well as the data.

library(coseq)
library(Biobase)
library(corrplot)

data("fietz")
counts <- exprs(fietz)
conds <- pData(fietz)$tissue ## Equivalent to the following: ## counts <- read.table("http://www.math.univ-toulouse.fr/~maugis/coseq/Fietz_mouse_counts.txt", ## header=TRUE) ## conds <- c(rep("CP",5),rep("SVZ",5),rep("VZ",5)) The coseq package fits either a Gaussian mixture model (Rau and Maugis-Rabusseau, 2017) or uses the K-means algorithm (Godichon-Baggioni et al., 2017) to cluster transformed normalized expression profiles of RNA-seq data. Normalized expression profiles correspond to the proportion of normalized reads observed for gene i with respect to the total observed for gene i across all samples: $p_{ij} = \frac{y_{ij}/s_{j} +1}{\sum_{j'} y_{ij'}/s_{j'} +1},$ where $$s_j$$ are normalization scaling factors (e.g., after applying TMM normalization to library sizes) and $$y_{ij}$$ represents the raw count for gene $$i$$ in sample $$j$$. ### Transformations for normalized profiles with the Gaussian mixture model Since the coordinates of $$\mathbf{p}_i$$ are linearly dependent (causing estimation problems for a Gaussian mixture distribution), weconsider either the arcsine or logit transformation of the normalized profiles $$p_{ij}$$: $g_{\text{arcsin}}(p_{ij}) = \text{arcsin}\left( \sqrt{ p_{ij} } \right) \in [0, \pi/2], \text{ and}$ $g_{\text{logit}}(p_{ij}) \text{log}_2 \left( \frac{p_{ij}}{1-p_{ij}} \right) \in (-\infty, \infty).$ Then the distribution of the transformed normalized expression profiles is modeled by a general multidimensional Gaussian mixture $f(.|\theta_K) = \sum_{k=1}^K \pi_k \phi(.|\mu_k,\Sigma_k)$ where $$\theta_K=(\pi_1,\ldots,\pi_{K-1},\mu_1,\ldots,\mu_K,\Sigma_1,\ldots,\Sigma_K)$$, $$\pi=(\pi_1,\ldots,\pi_K)$$ are the mixing proportions and $$\phi(.|\mu_k,\Sigma_k)$$ is the $$q$$-dimensional Gaussian density function with mean $$\mu_k$$ and covariance matrix $$\Sigma_k$$. To estimate mixture parameters $$\theta_K$$ by computing the maximum likelihood estimate (MLE), an Expectation-Maximization (EM) algorithm is used via the Rmixmod package. Finally, let $$\hat t_{ik}$$ be the conditional probability that observation $$i$$ arises from the $$k$$th component of the mixture $$f(.|\hat \theta_{\hat K})$$. Each observation $$i$$ is assigned to the component maximizing the conditional probability $$\hat t_{ik}$$ i.e., using the so-called maximum a posteriori (MAP) rule. ### Transformations for normalized profiles with K-means For the K-means algorithm, we consider three separate transformations of the profiles $$p_{i}$$ that are well adapted to compositional data (see Godichon-Baggioni et al., 2017 for more details): • Identity (i.e., no transformation) • Centered log ratio (CLR) transformation • Log centered log ratio (logCLR) transformation Then, the aim is to minimize the sum of squared errors (SSE), defined for each set of clustering $$\mathcal{C}^{(K)}=\left( C_{1},...,C_{k}\right)$$ by $\begin{equation*} \text{SSE} \left( \mathcal{C}^{(K)}\right) := \sum_{k=1}^{K}\sum_{i \in C_{k}} \left\| h \left( y_{i}\right) - \mu_{k,h} \right\|_{2}^{2} , \end{equation*}$ with $$i \in C_{k}$$ if $$\left\| h\left( y_{i}\right) - \mu_{k,h} \right\|_{2} = \min_{k'=1,\ldots,K} \left\| y_{i}- \mu_{k',h} \right\|_{2}$$, and $\mu_{k,h}= \frac{1}{|C_{k} |}\sum_{i \in C_{k}}h \left( y_{i} \right),$ and $$h$$ is the chosen transformation. In order to minimize the SSE, we use the well-known MacQueen’s K-means algorithm. ### Model selection Because the number of clusters $$K$$ is not known a priori, we fit a collection of models (here $$K$$ = 2, …, 25) and use either the Integrated Completed Likelihood (ICL) criterion (in the case of the Gaussian mixture model) or the slope heuristics (in the case of the K-means algorithm) to select the best model in terms of fit, complexity, and cluster separation. ### Other options If desired, we can set a filtering cutoff on the mean normalized counts via the meanFilterCutoff argument to remove very weakly expressed genes from the co-expression analysis; in the interest of faster computation for this vignette, in this example we filter all genes with mean normalized expression less than 200 (for the Gaussian mixture model) or less than 50 (for the K-means algorithm). Note that if desired, parallel execution using BiocParallel can be specified via the parallel argument: run <- coseq(..., parallel=TRUE) ## Running coseq The collection of co-expression models for the logCLR-transformed profiles using the K-means algorithm may be obtained as follows (note that we artificially set the number of maximum allowed iterations and number of random starts to be low for faster computational time in this vignette): runLogCLR <- coseq(counts, K=2:25, transformation="logclr",norm="TMM", meanFilterCutoff=50, model="kmeans", nstart=1, iter.max=10) ## **************************************** ## coseq analysis: kmeans approach & logclr transformation ## K = 2 to 25 ## **************************************** ## Running K = 2 ... ## Running K = 3 ... ## Running K = 4 ... ## Running K = 5 ... ## Running K = 6 ... ## Running K = 7 ... ## Running K = 8 ... ## Running K = 9 ... ## Running K = 10 ... ## Running K = 11 ... ## Running K = 12 ... ## Running K = 13 ... ## Running K = 14 ... ## Running K = 15 ... ## Running K = 16 ... ## Running K = 17 ... ## Running K = 18 ... ## Running K = 19 ... ## Running K = 20 ... ## Running K = 21 ... ## Running K = 22 ... ## Running K = 23 ... ## Running K = 24 ... ## Running K = 25 ... The collection of Gaussian mixture models for the arcsine-transformed and logit-transformed profiles may be obtained as follows (as before, we set the number of iterations to be quite low for computational speed in this vignette): runArcsin <- coseq(counts, K=2:20, model="Normal", transformation="arcsin", meanFilterCutoff=200, iter=10) ## **************************************** ## coseq analysis: Normal approach & arcsin transformation ## K = 2 to 20 ## **************************************** ## Running K = 2 ... ## Running K = 3 ... ## Running K = 4 ... ## Running K = 5 ... ## Running K = 6 ... ## Running K = 7 ... ## Running K = 8 ... ## Running K = 9 ... ## Running K = 10 ... ## Running K = 11 ... ## Running K = 12 ... ## Running K = 13 ... ## Running K = 14 ... ## Running K = 15 ... ## Running K = 16 ... ## Running K = 17 ... ## Running K = 18 ... ## Running K = 19 ... ## Running K = 20 ... runLogit <- coseq(counts, K=2:20, model="Normal", transformation="logit", meanFilterCutoff=200, verbose=FALSE, iter=10) ## **************************************** ## coseq analysis: Normal approach & logit transformation ## K = 2 to 20 ## **************************************** In all cases, the resulting output of a call to coseq is an S4 object of class coseqResults. class(runArcsin) ## [1] "coseqResults" ## attr(,"package") ## [1] "coseq" runArcsin ## An object of class coseqResults ## 4230 features by 15 samples. ## Models fit: K = 2 ... 20 ## Chosen clustering model: K = 8 To choose the most appropriate transformation to use (arcsine versus logit) in a Gaussian mixture model, we may use the corrected ICL, where the minimum value corresponds to the selected model. Note that this feature is not available for the K-means algorithm. compareICL(list(runArcsin, runLogit)) This indicates that the preferred transformation is the arcsine, and the selected model has$K = ${r} ncol(runArcsin) clusters. We can additionally explore how similar the models with$K = ${r ncol(runArcsin)-1} to {r ncol(runArcsin)+1} are using the adjusted Rand index (ARI) via the compareARI function. Values close to 1 indicate perfect agreement, while values close to 0 indicate near random partitions. compareARI(runArcsin, K=8:12) ## K=8 K=9 K=10 K=11 K=12 ## K=8 1 0.36 0.5 0.43 0.51 ## K=9 1 0.44 0.4 0.46 ## K=10 1 0.39 0.39 ## K=11 1 0.38 ## K=12 1 Note that because the results of coseq depend on the initialization point, results from one run to another may vary; as such, in practice, it is typically a good idea to re-estimate the same collection of models a few times (e.g., 5) to avoid problems linked to initialization. ## Exploring coseq results Results from a coseq analysis can be explored and summarized in a variety of ways. First, a call to the summary function provides the number of clusters selected for the ICL model selection approach, number of genes assigned to each cluster, and if desired the per-gene cluster means. summary(runArcsin) ## ************************************************* ## Model: Gaussian_pk_Lk_Ck ## Transformation: arcsin ## ************************************************* ## Clusters fit: 2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20 ## Clusters with errors: --- ## Selected number of clusters via ICL: 8 ## ICL of selected model: -365687.2 ## ************************************************* ## Number of clusters = 8 ## ICL = -365687.2 ## ************************************************* ## Cluster sizes: ## Cluster 1 Cluster 2 Cluster 3 Cluster 4 Cluster 5 Cluster 6 Cluster 7 ## 338 1047 259 418 945 357 177 ## Cluster 8 ## 689 ## ## Number of observations with MAP > 0.90 (% of total): ## 3429 (81.06%) ## ## Number of observations with MAP > 0.90 per cluster (% of total per cluster): ## Cluster 1 Cluster 2 Cluster 3 Cluster 4 Cluster 5 Cluster 6 Cluster 7 ## 252 809 249 397 713 284 155 ## (74.56%) (77.27%) (96.14%) (94.98%) (75.45%) (79.55%) (87.57%) ## Cluster 8 ## 570 ## (82.73%) Next, a variety of plots may be explored using the plot function: • logLike (log-likelihood plotted versus number of clusters), • ICL (ICL plotted versus number of clusters), • profiles(line plots of profiles in each cluster), • boxplots (boxplots of profiles in each cluster), • probapost_boxplots (boxplots of maximum conditional probabilities per cluster), • probapost_barplots (number of observations with a maximum conditional probability greater than a given threshold per cluster), • probapost_histogram (histogram of maximum conditional probabilities over all clusters). By default, all of these plots are produced simultaneously unless specific graphics are requested via the graphs argument. In addition, a variety of options are available to specify the number of graphs per row/column in a grid layout, whether profiles should be averaged over conditions, whether condition labels should be used, whether clusters should be ordered according to their similarity, etc. Note that the histogram of maximum conditional probabilities of cluster membership for all genes (probapost_histogram), and boxplots and barplots of maximum conditional probabilities of cluster membership for the genes assigned to each cluster (probapost_boxplots, probapost_barplots) help to evaluate the degree of certitude accorded by the model in assigning genes to clusters, as well as whether some clusters are attribued a greater degree of uncertainty than others. ## To obtain all plots ## plot(runArcsin) plot(runArcsin, graphs="boxplots") ##$boxplots

plot(runArcsin, graphs="boxplots", conds=conds)
## $boxplots plot(runArcsin, graphs="boxplots", conds=conds, average_over_conds=TRUE) ##$boxplots

plot(runArcsin, graphs="profiles")
## $profiles plot(runArcsin, graphs="probapost_boxplots") ##$probapost_boxplots

plot(runArcsin, graphs="probapost_histogram")
## $probapost_histogram If desired the per-cluster correlation matrices estimated in the Gaussian mixture model may be obtained by calling NormMixParam. These matrices may easily be visualized using the corrplot package. rho <- NormMixParam(runArcsin)$rho
## Covariance matrix for cluster 1
rho1 <- rho[,,1]
colnames(rho1) <- rownames(rho1) <- paste0(colnames(rho1), "\n", conds)
corrplot(rho1)

Finally, cluster labels (as well as a variety of other information) are easily accessible via a set of accessor functions.

labels <- clusters(runArcsin)
table(labels)
## labels
##    1    2    3    4    5    6    7    8
##  338 1047  259  418  945  357  177  689
metadata(runArcsin)
## $nbCluster ## K=2 K=3 K=4 K=5 K=6 K=7 K=8 K=9 K=10 K=11 K=12 K=13 K=14 K=15 K=16 ## 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 ## K=17 K=18 K=19 K=20 ## 17 18 19 20 ## ##$logLike
##      K=2      K=3      K=4      K=5      K=6      K=7      K=8      K=9
## 180272.9 182045.0 183805.2 185539.8 186043.9 187322.1 188077.6 188570.5
##     K=10     K=11     K=12     K=13     K=14     K=15     K=16     K=17
## 189187.2 189640.5 189932.5 190433.9 190602.6 190838.9 191366.7 191856.1
##     K=18     K=19     K=20
## 192232.4 192114.0 192388.9
##
## $ICL ## K=2 K=3 K=4 K=5 K=6 K=7 K=8 ## -357870.3 -360214.1 -362198.8 -364266.4 -364421.6 -365116.5 -365687.2 ## K=9 K=10 K=11 K=12 K=13 K=14 K=15 ## -365301.7 -365167.8 -365052.6 -364535.1 -364489.9 -363622.3 -363142.9 ## K=16 K=17 K=18 K=19 K=20 ## -362662.5 -362427.2 -362234.1 -360745.0 -360136.2 ## ##$nbClusterError
## integer(0)
##