Email updates

Keep up to date with the latest news and content from BMC Bioinformatics and BioMed Central.

Open Access Highly Accessed Research article

Unsupervised Bayesian linear unmixing of gene expression microarrays

Cécile Bazot1*, Nicolas Dobigeon1*, Jean-Yves Tourneret1, Aimee K Zaas2, Geoffrey S Ginsburg2 and Alfred O Hero III3

Author Affiliations

1 University of Toulouse, IRIT/INP-ENSEEIHT, 2 rue Camichel, 31071 Toulouse cedex 7, BP 7122, France

2 Department of Medicine and Institute for Genome Sciences and Policy, Duke University, Durham, North Carolina, USA

3 Center for Computational Biology and Bioinformatics and EECS Department, University of Michigan, 1301 Beal Avenue, Ann Arbor, MI, 48109-2122, USA

For all author emails, please log on.

BMC Bioinformatics 2013, 14:99  doi:10.1186/1471-2105-14-99


The electronic version of this article is the complete one and can be found online at: http://www.biomedcentral.com/1471-2105/14/99


Received:5 December 2012
Accepted:24 January 2013
Published:19 March 2013

© 2013 Bazot et al.; licensee BioMed Central Ltd.

This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract

Background

This paper introduces a new constrained model and the corresponding algorithm, called unsupervised Bayesian linear unmixing (uBLU), to identify biological signatures from high dimensional assays like gene expression microarrays. The basis for uBLU is a Bayesian model for the data samples which are represented as an additive mixture of random positive gene signatures, called factors, with random positive mixing coefficients, called factor scores, that specify the relative contribution of each signature to a specific sample. The particularity of the proposed method is that uBLU constrains the factor loadings to be non-negative and the factor scores to be probability distributions over the factors. Furthermore, it also provides estimates of the number of factors. A Gibbs sampling strategy is adopted here to generate random samples according to the posterior distribution of the factors, factor scores, and number of factors. These samples are then used to estimate all the unknown parameters.

Results

Firstly, the proposed uBLU method is applied to several simulated datasets with known ground truth and compared with previous factor decomposition methods, such as principal component analysis (PCA), non negative matrix factorization (NMF), Bayesian factor regression modeling (BFRM), and the gradient-based algorithm for general matrix factorization (GB-GMF). Secondly, we illustrate the application of uBLU on a real time-evolving gene expression dataset from a recent viral challenge study in which individuals have been inoculated with influenza A/H3N2/Wisconsin. We show that the uBLU method significantly outperforms the other methods on the simulated and real data sets considered here.

Conclusions

The results obtained on synthetic and real data illustrate the accuracy of the proposed uBLU method when compared to other factor decomposition methods from the literature (PCA, NMF, BFRM, and GB-GMF). The uBLU method identifies an inflammatory component closely associated with clinical symptom scores collected during the study. Using a constrained model allows recovery of all the inflammatory genes in a single factor.

Background

Factor analysis methods such as principal component analysis (PCA) have been widely studied and can be used for discovering the patterns of differential expression in time course and/or multiple treatment biological experiments using gene or protein microarray samples. These methods aim at finding a decomposition of the observation matrix <a onClick="popup('http://www.biomedcentral.com/1471-2105/14/99/mathml/M1','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/14/99/mathml/M1">View MathML</a> whose rows (respectively columns) are indexed by gene index (respectively sample index). Typically, in gene expression analysis, the number N of samples is much less than the number G of genes. For example, in an Affymetrix HU133 gene chip, the number G of genes may range from ten to twenty thousand depending on the type of chip description file (CDF) processing used to translate the oligonucleotide fragments to gene labels whereas we only analyze about a hundred of samples.

This decomposition expresses each of the N samples as a particular linear combination of R characteristic gene expression signatures, also referred to as factors, with appropriate proportions (or contributions), called factor scores, following a linear mixing model

<a onClick="popup('http://www.biomedcentral.com/1471-2105/14/99/mathml/M2','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/14/99/mathml/M2">View MathML</a>

(1)

where <a onClick="popup('http://www.biomedcentral.com/1471-2105/14/99/mathml/M3','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/14/99/mathml/M3">View MathML</a> represents the factor loading matrix, <a onClick="popup('http://www.biomedcentral.com/1471-2105/14/99/mathml/M4','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/14/99/mathml/M4">View MathML</a> the factor score matrix and <a onClick="popup('http://www.biomedcentral.com/1471-2105/14/99/mathml/M5','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/14/99/mathml/M5">View MathML</a> is a matrix containing noise samples. Each sample yi (i = 1, …, N), corresponding to the i-th column of the observed gene expression matrix Y, is a vector of G gene expression levels that can be expressed as

<a onClick="popup('http://www.biomedcentral.com/1471-2105/14/99/mathml/M6','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/14/99/mathml/M6">View MathML</a>

(2)

where mr is the r-th column of M, ar, i denotes the (r, i)-th element of the matrix A, ai and ni are the i-th column of A and N respectively. The number of factors R in the decomposition is usually much less than the number of samples N. Traditional factor analysis methods such as PCA require the experimenter to specify the desired number of factors to be estimated. However, some recent Bayesian factor analysis methods are totally unsupervised in the sense that the number of factors is directly estimated from the data [1-3].

The model (1) is identical to the standard factor analysis model [4] for which the columns of M are called factors and should correspond to biological signatures (or pathways). Note that the elements of the matrix M are referred to as factor loadings, and the columns of A are the factor scores. Approaches to fitting the model (1) to data include principal component analysis [5,6], least squares matrix factorization [7,8], finite mixture modeling [9,10], and Bayesian factor analysis [4,11,12].

This paper presents a new Bayesian factor analysis method called unsupervised Bayesian linear unmixing (uBLU), that estimates the number of factors and incorporates non-negativity constraints on the factors and factor scores, as well as a sum-to-one constraint for the factor scores. The uBLU method presented here differs from the BLU method, developed in [13] for hyperspectral imaging and applied to gene microarray expression analysis in [14]. Note that BLU requires user specification of the number of factors while uBLU determines the number of factors using Bayesian birth-death model. The positivity and sum-to-one constraints are natural in gene microarray analysis when the entries of the observation matrix are non-negative and when a proportional decomposition is desired. Thus each factor score corresponds to the concentration (or proportion) of a particular factor to a given sample. The advantage of this representation for gene expression analysis is twofold: i) the factor scores correspond to the strengths of each gene contributing to that factor; ii) for each gene chip the factor scores give the relative abundance of each factor present in the chip. For example, a gene having a large loading level (close to one) for a particular factor should have a small loading (close to zero) for all other factors. In this way, as opposed to other factor analysis methods, there is less multiplexing making it easier to associate specific genes to specific factors and vice versa.

A similar approach, based on NMR spectral imaging and called the Bayesian decomposition (BD), has been previously developed by Moloshok et al. and applied to gene expression data [11]. More recently, the coordinated gene activity in pattern sets method (CoGAPS), available as an open R-source [12], combines the GAPS-MCMC matrix factorization algorithm with a threshold-independent statistic to infer activity in specific gene sets. However, these approaches require cold restarts of the algorithm with different number of factors and with different random seeds to avoid the large number of local minima encountered in the process of fitting the matrix factorization model MA to the data Y. In contrast, the proposed uBLU algorithm uses a judicious model to reduce sensitivity to local minima rather than using cold restarts. The novelty of the uBLU model is that it consists of: (1) a birth-death process to infer the number of factors; (2) a positivity constraint on the loading and score matrices M, A to restrict the space of solutions; (3) a sum-to-one constraint on the columns of A to further restrict the solution space. The uBLU model is justified for non-negative data problems like gene expression analysis and produces an estimate of the non-negative factors in addition to their proportional representation in each sample.

Bayesian linear unmixing, traditionally used for hyperspectral image analysis (see [13] for example), is one of many possible factor analysis methods that could be applied to gene expression analysis. These methods include non-negative matrix factorization (NMF) [7,8], independent component analysis (ICA) [15], Bayesian decomposition (BD) [11], PCA [5], bi-clustering [16], penalized matrix decomposition (PMD) [2], Bayesian factor regression modeling (BFRM) [1], and more recently the gradient-based algorithm of Nikulin et al. for general matrix factorization (GB-GMF) [17]. Contrary to uBLU, the PCA, ICA, BFRM, GB-GMF, bi-clustering and PMD methods do not account for non-negativity of the factor loadings and factor scores. On the other hand, NMF does not account for sum-to-one constraints on the columns of the factor score matrix. Contrary to PCA and ICA, uBLU does not impose orthogonality or independence on the factors, as well as the GB-GMF algorithm. These relaxed assumptions might better represent what is known about the preponderance of overlap and dependency in biological pathways. Finally, uBLU naturally accommodates Bayesian estimation of the number of factors, like BFRM. Note that BFRM has been specifically developed for gene expression data [1].

In this paper we provide comparative studies that establish quantitative performance advantages of the proposed constrained model and its corresponding uBLU algorithm with respect to PCA, NMF, BFRM and GB-GMF for time-varying gene expression analysis, using synthetic data with known ground truth. We also illustrate the application of uBLU to the analysis of a real gene expression dataset from a recent viral challenge study [18] in which several subjects were administered viral inoculum and gene expression time course data were collected over a period of several days. Using these data, we may infer relationships between genes and symptoms and examine how the human host response to viral infection evolves with time.

Methods

Mathematical constrained model

Let yi represent a gene microarray vector of G gene expression levels. The elements of yi have units of hybridization abundance levels with non-negative values. In the context of gene expression data, the starting point for Bayesian linear unmixing is the linear mixing model (LMM)

<a onClick="popup('http://www.biomedcentral.com/1471-2105/14/99/mathml/M7','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/14/99/mathml/M7">View MathML</a>

(3)

where R is the number of distinct gene signatures that can be present in the chip, mr = [m1, r, …, mG, r]T is the r-th gene signature vector, mg, r ≥ 0 is the strength of the g-th gene (g = 1, …, G) in the r-th signature (r = 1, …, R), and ar, i is the relative contribution of the r-th signature vector to the i-th sample yi, where ar,i ∈ [0, 1] and <a onClick="popup('http://www.biomedcentral.com/1471-2105/14/99/mathml/M8','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/14/99/mathml/M8">View MathML</a>. Here ni denotes the residual error of the LMM representation. For a matrix of N data samples <a onClick="popup('http://www.biomedcentral.com/1471-2105/14/99/mathml/M9','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/14/99/mathml/M9">View MathML</a>, the LMM can be rewritten with matrix notations

<a onClick="popup('http://www.biomedcentral.com/1471-2105/14/99/mathml/M10','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/14/99/mathml/M10">View MathML</a>

(4)

where <a onClick="popup('http://www.biomedcentral.com/1471-2105/14/99/mathml/M11','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/14/99/mathml/M11">View MathML</a>, <a onClick="popup('http://www.biomedcentral.com/1471-2105/14/99/mathml/M12','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/14/99/mathml/M12">View MathML</a> and <a onClick="popup('http://www.biomedcentral.com/1471-2105/14/99/mathml/M13','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/14/99/mathml/M13">View MathML</a> represent the factor score matrix, the factor loading matrix and the noise matrix, respectively. The matrices M, A satisfy positivity and sum-to-one constraints defined by

<a onClick="popup('http://www.biomedcentral.com/1471-2105/14/99/mathml/M14','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/14/99/mathml/M14">View MathML</a>

(5)

where mg, r denotes the (g, r)-th element of matrix M. The constraints (5) arise naturally when dealing with positive data for which one is seeking the relative contribution of positive factors that have the same numerical characteristics as the data, i.e., the signature mr is itself interpretable as a vector of hybridization abundances.

The objective of linear unmixing is to simultaneously estimate the factor matrix M and the factor score matrix A from the available N data samples. The representation (1) is rank deficient since A has rank N − 1. This presents algorithmic challenges for solving the unmixing problem. Several algorithms have been proposed in the context of hyperspectral imaging to solve similar problems [6,19]. Most of these algorithms perform unmixing in a two step procedure where M is estimated first using an endmember extraction algorithm (EEA) followed by a constrained linear least squares step to estimate A. A common (but restrictive) assumption in these algorithms is that some samples in the dataset are “pure” in the sense that the linear combination of (2) contains a unique factor, say mr, with factor score ar, i. Recently, this assumption has been relaxed by applying a hierarchical Bayesian approach, called Bayesian linear unmixing (BLU). The resulting algorithm requires the number R of factors to be specified (see [13] for details). Here we extend BLU to a fully unsupervised algorithm, called unsupervised BLU (uBLU), that estimates R using a birth-death model and a Gibbs sampler. The Gibbs sampler produces an estimate of the entire joint posterior distribution of the model parameters, resulting in a fully Bayesian estimator of the number of factors R, the factor loadings M, and the factor scores A. The uBLU model is described in the next subsection and the Gibbs sampling algorithm is given in the Appendix. In the Results and discussion section below we demonstrate the performance advantages of uBLU as a factor analysis model for simulated and real gene expression data.

Unsupervised Bayesian linear unmixing algorithm

The BLU algorithm studied in [13] generates samples distributed according to the posterior distribution of M and A given the number R of factors for appropriate prior distributions assigned to the mixing parameters in (2). First, the residual errors ni in (2) are assumed to be independent identically distributed (i.i.d.) according to zero-mean Gaussian distributions: <a onClick="popup('http://www.biomedcentral.com/1471-2105/14/99/mathml/M15','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/14/99/mathml/M15">View MathML</a> for i = 1, …, N, where IG denotes the identity matrix of dimension G × G.

The number of factors R to be estimated by the proposed uBLU algorithm is assigned a discrete uniform prior distribution on [2, …, Rmax]

<a onClick="popup('http://www.biomedcentral.com/1471-2105/14/99/mathml/M16','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/14/99/mathml/M16">View MathML</a>

(6)

where Rmax is the maximal number of factors present in the mixture.

Because of the constraints in (5), the data samples yi (i = 1, …, N) live in a lower-dimensional subspace of <a onClick="popup('http://www.biomedcentral.com/1471-2105/14/99/mathml/M17','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/14/99/mathml/M17">View MathML</a> (whose dimension is upper-bounded by K − 1) denoted as <a onClick="popup('http://www.biomedcentral.com/1471-2105/14/99/mathml/M18','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/14/99/mathml/M18">View MathML</a> (Rmax − 1 ≤ K ≤ G). This subspace can be identified by a standard dimension reduction procedure, such as PCA. Hence, instead of estimating the factor loadings <a onClick="popup('http://www.biomedcentral.com/1471-2105/14/99/mathml/M19','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/14/99/mathml/M19">View MathML</a> (r = 1, …, R), we propose to estimate their corresponding projections <a onClick="popup('http://www.biomedcentral.com/1471-2105/14/99/mathml/M20','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/14/99/mathml/M20">View MathML</a> onto this subspace. Specifically, these projections can be represented as

<a onClick="popup('http://www.biomedcentral.com/1471-2105/14/99/mathml/M21','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/14/99/mathml/M21">View MathML</a>

(7)

where <a onClick="popup('http://www.biomedcentral.com/1471-2105/14/99/mathml/M22','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/14/99/mathml/M22">View MathML</a> is the empirical mean of the data matrix Y and P is the (K − 1) × G appropriate projection matrix that projects onto <a onClick="popup('http://www.biomedcentral.com/1471-2105/14/99/mathml/M23','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/14/99/mathml/M23">View MathML</a>, which can be constructed from the principal eigenvectors of the empirical covariance matrix of Y. This dimension reduction procedure allows us to work in a lower-dimensional subspace without loss of information, and reduces significantly the computational complexity of the Bayes estimator of the factor loadings. A multivariate Gaussian distribution (MGD) truncated on a subset <a onClick="popup('http://www.biomedcentral.com/1471-2105/14/99/mathml/M24','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/14/99/mathml/M24">View MathML</a> is chosen as prior distribution for the projected factors tr. The subset <a onClick="popup('http://www.biomedcentral.com/1471-2105/14/99/mathml/M25','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/14/99/mathml/M25">View MathML</a> is defined in order to ensure the non-negativity constraint on mr (see [13])

<a onClick="popup('http://www.biomedcentral.com/1471-2105/14/99/mathml/M26','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/14/99/mathml/M26">View MathML</a>

(8)

More precisely, <a onClick="popup('http://www.biomedcentral.com/1471-2105/14/99/mathml/M27','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/14/99/mathml/M27">View MathML</a> is obtained by noting that <a onClick="popup('http://www.biomedcentral.com/1471-2105/14/99/mathml/M28','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/14/99/mathml/M28">View MathML</a> and by looking for the vectors tr such that all components of <a onClick="popup('http://www.biomedcentral.com/1471-2105/14/99/mathml/M29','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/14/99/mathml/M29">View MathML</a> are non-negative. To estimate the mean vectors er of these truncated MGDs, one can use a standard endmember extraction algorithm (EEA) common in hyperspectral imaging, e.g. N-FINDR [19]. To summarize, the prior distribution for the projected factor tr is

<a onClick="popup('http://www.biomedcentral.com/1471-2105/14/99/mathml/M30','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/14/99/mathml/M30">View MathML</a>

(9)

where <a onClick="popup('http://www.biomedcentral.com/1471-2105/14/99/mathml/M31','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/14/99/mathml/M31">View MathML</a> denotes the truncated MGD with mean vector er and covariance matrix <a onClick="popup('http://www.biomedcentral.com/1471-2105/14/99/mathml/M32','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/14/99/mathml/M32">View MathML</a>, with <a onClick="popup('http://www.biomedcentral.com/1471-2105/14/99/mathml/M33','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/14/99/mathml/M33">View MathML</a> a fixed hyperparameter. Assuming the vectors tr, for r = 1, …, R, are a priori independent, the prior distribution for the projected factor matrix T = [t1,…,tR] is

<a onClick="popup('http://www.biomedcentral.com/1471-2105/14/99/mathml/M34','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/14/99/mathml/M34">View MathML</a>

(10)

where ∝ stands for “proportional to”, ∥·∥ is the standard l2-norm, <a onClick="popup('http://www.biomedcentral.com/1471-2105/14/99/mathml/M35','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/14/99/mathml/M35">View MathML</a> denotes the indicator function on the set <a onClick="popup('http://www.biomedcentral.com/1471-2105/14/99/mathml/M36','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/14/99/mathml/M36">View MathML</a>, E = [e1, …, eR] and <a onClick="popup('http://www.biomedcentral.com/1471-2105/14/99/mathml/M37','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/14/99/mathml/M37">View MathML</a>.

The sum-to-one constraint for the factor scores ai, for each observed sample i (i = 1, …, N), allows this vector ai to be rewritten as

<a onClick="popup('http://www.biomedcentral.com/1471-2105/14/99/mathml/M38','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/14/99/mathml/M38">View MathML</a>

(11)

and <a onClick="popup('http://www.biomedcentral.com/1471-2105/14/99/mathml/M39','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/14/99/mathml/M39">View MathML</a>. Note here that any component of ai could be expressed as a function of the others, i.e., <a onClick="popup('http://www.biomedcentral.com/1471-2105/14/99/mathml/M40','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/14/99/mathml/M40">View MathML</a>. The last component aR, i has been chosen here for notation simplicity. To ensure the positivity constraint, the subvectors a1:R − 1, i must belong to the simplex

<a onClick="popup('http://www.biomedcentral.com/1471-2105/14/99/mathml/M41','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/14/99/mathml/M41">View MathML</a>

(12)

where ∥·∥1 is the l1 norm (<a onClick="popup('http://www.biomedcentral.com/1471-2105/14/99/mathml/M42','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/14/99/mathml/M42">View MathML</a>) and ai0 stands for the set of inequalities {ar,i ≥ 0}r = 1, …, R. Following the model in [13], we propose to assign uniform distributions over the simplex <a onClick="popup('http://www.biomedcentral.com/1471-2105/14/99/mathml/M43','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/14/99/mathml/M43">View MathML</a> as priors for the subvectors a1:R − 1, i (i = 1, …, N), i.e.,

<a onClick="popup('http://www.biomedcentral.com/1471-2105/14/99/mathml/M44','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/14/99/mathml/M44">View MathML</a>

(13)

For the prior distribution on the variance σ2 of the residual errors we chose a conjugate inverse-Gamma distribution with parameters ν / 2 and γ / 2

<a onClick="popup('http://www.biomedcentral.com/1471-2105/14/99/mathml/M45','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/14/99/mathml/M45">View MathML</a>

(14)

The shape parameter ν is a fixed hyperparameter whereas the scale parameter γ will be adjustable, as in [13]. A non-informative Jeffreys’ prior is chosen as prior distribution for the hyperparameter γ, i.e.,

<a onClick="popup('http://www.biomedcentral.com/1471-2105/14/99/mathml/M46','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/14/99/mathml/M46">View MathML</a>

(15)

The resulting hierarchical structure of the proposed uBLU model is summarized in the directed acyclic graph (DAG) presented in Additional file 1: Figure S1.

Additional file 1. Supplementary materials on algorithm details and performance validation. Directed acyclic graph (DAG) of the model and flowchart of the proposed algorithm are provided in this additional file. More results on synthetic datasets are also presented to validate the proposed Bayesian algorithm, including a convergence diagnosis.

Format: PDF Size: 774KB Download file

This file can be viewed with: Adobe Acrobat ReaderOpen Data

The model defined in (1) and the Gaussian assumption for the noise vectors n1, …, nN allow the likelihood of y1, …, yN to be determined

<a onClick="popup('http://www.biomedcentral.com/1471-2105/14/99/mathml/M47','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/14/99/mathml/M47">View MathML</a>

(16)

Multiplying this likelihood by the parameter priors defined in (10), (13), (14) and (6), and integrating out the nuisance parameter γ, the posterior distribution of the unknown parameter vector Θ = {M, A, σ2, R} can be expressed as

<a onClick="popup('http://www.biomedcentral.com/1471-2105/14/99/mathml/M48','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/14/99/mathml/M48">View MathML</a>

(17)

Considering the parameters to be a priori independent, the following result can be obtained

<a onClick="popup('http://www.biomedcentral.com/1471-2105/14/99/mathml/M49','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/14/99/mathml/M49">View MathML</a>

(18)

where f(A|R), f(T|E, s2, R) and f(σ2|ν, γ) are respectively the prior distributions of the factor score matrix A, the projected factor matrix T and the noise variance σ2 previously defined.

Due to the constraints enforced on the data, the posterior distribution f(M, A, R|Y) obtained from the proposed hierarchical structure is too complex to derive analytical expressions of the Bayesian estimators, e.g., the minimum mean square (MMSE) and maximum a posteriori (MAP) estimators. In such case, it is natural to use Markov chain Monte Carlo (MCMC) methods [20] to generate samples M(t), A(t) and R(t) asymptotically distributed according to f(M, A, R|Y). However, the dimensions of the factor loading matrix M and the factor score matrix A depend on the unknown number R of signatures to be identified. As a consequence, sampling from f(M, A, R|Y) requires exploring parameter spaces of different dimensions. To solve this dimension matching problem, we include a birth/death process within the MCMC procedure. Specifically, a birth, death or switch move is chosen at each iteration of the algorithm (see the Appendix and [21]). This birth-death model differs from the classical reversible-jump MCMC (RJ-MCMC) (as defined in [21]) in the sense that, for the birth-death model, each move is accepted or rejected at each iteration using the likelihood ratio between the current state and the new state proposed by the algorithm. The factor matrix M, the factor score matrix A and the noise variance σ2 are then updated, conditionally upon the number of factors R, using Gibbs moves.

After a sufficient number of iterations (Nmc iterations, including a burn-in period of Nbi iterations), the traditional Bayesian estimators (e.g., MMSE and MAP) can be approximated using the generated samples M(t), A(t) and R(t). First, the generated samples are used to approximate the MAP estimator of the number of factors

<a onClick="popup('http://www.biomedcentral.com/1471-2105/14/99/mathml/M50','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/14/99/mathml/M50">View MathML</a>

(19)

where Nk is the number of generated samples <a onClick="popup('http://www.biomedcentral.com/1471-2105/14/99/mathml/M51','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/14/99/mathml/M51">View MathML</a> satisfying R(t) = k and <a onClick="popup('http://www.biomedcentral.com/1471-2105/14/99/mathml/M52','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/14/99/mathml/M52">View MathML</a>. Then, conditioned on <a onClick="popup('http://www.biomedcentral.com/1471-2105/14/99/mathml/M53','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/14/99/mathml/M53">View MathML</a>, the joint MAP estimator <a onClick="popup('http://www.biomedcentral.com/1471-2105/14/99/mathml/M54','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/14/99/mathml/M54">View MathML</a> of the factor and factor score matrices is determined as follows

<a onClick="popup('http://www.biomedcentral.com/1471-2105/14/99/mathml/M55','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/14/99/mathml/M55">View MathML</a>

(20)

Results and discussion

The proposed method consists of estimating simultaneously the matrices M and A defined in (1), under the positivity and sum-to-one constraints mentioned previously, in a fully unsupervised framework, i.e., the number of factors R is also estimated from the data. A Gibbs sampler algorithm is designed that generates samples distributed according to the posterior distribution associated to the proposed uBLU model. For more details about the Gibbs sampling strategy, see the Appendix.

Simulations on synthetic data

To illustrate the performance of the proposed Bayesian factor decomposition, we first present simulations conducted on synthetic data. More extensive simulation results are reported in the Additional file 1.

Simulation scenario

Several synthetic datasets <a onClick="popup('http://www.biomedcentral.com/1471-2105/14/99/mathml/M56','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/14/99/mathml/M56">View MathML</a> were generated. The experiments presented here correspond to the expression values of G = 512 genes (for datasets <a onClick="popup('http://www.biomedcentral.com/1471-2105/14/99/mathml/M57','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/14/99/mathml/M57">View MathML</a>, <a onClick="popup('http://www.biomedcentral.com/1471-2105/14/99/mathml/M58','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/14/99/mathml/M58">View MathML</a> and <a onClick="popup('http://www.biomedcentral.com/1471-2105/14/99/mathml/M59','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/14/99/mathml/M59">View MathML</a>) or G=12000 genes (for dataset <a onClick="popup('http://www.biomedcentral.com/1471-2105/14/99/mathml/M60','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/14/99/mathml/M60">View MathML</a>) with N = 128 samples. Each sample is composed of exactly R=3 factors mixed using the linear mixing model in (1). The factors of the first dataset <a onClick="popup('http://www.biomedcentral.com/1471-2105/14/99/mathml/M61','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/14/99/mathml/M61">View MathML</a> have been generated so that only a few genes affect each factor. For the second dataset <a onClick="popup('http://www.biomedcentral.com/1471-2105/14/99/mathml/M62','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/14/99/mathml/M62">View MathML</a>, realistic factors have been extracted from real genetic datasets. The third dataset <a onClick="popup('http://www.biomedcentral.com/1471-2105/14/99/mathml/M63','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/14/99/mathml/M63">View MathML</a> has been generated enforcing the factors to be orthogonal but not necessarily positive whereas in the forth dataset, <a onClick="popup('http://www.biomedcentral.com/1471-2105/14/99/mathml/M64','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/14/99/mathml/M64">View MathML</a>, factors are orthogonal and positive. These simulation conditions are summarized in Table 1.

Table 1. Synthetic datasets<a onClick="popup('http://www.biomedcentral.com/1471-2105/14/99/mathml/M65','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/14/99/mathml/M65">View MathML</a>

In each case, the R = 3 factors were mixed in random proportions (factor scores), with positivity and sum-to-one constraints. All synthetic datasets were corrupted by an i.i.d. Gaussian noise sequence. The signal-to-noise ratio is SNRi = 20 dB where<a onClick="popup('http://www.biomedcentral.com/1471-2105/14/99/mathml/M66','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/14/99/mathml/M66">View MathML</a> for each sample i(i = 1, …, N).

Proposed method (uBLU)

The first step of the algorithm consists of estimating the number of factors R involved in the mixture, and hence determining the dimensions of the matrices M and A, using the maximum a posteriori (MAP) estimator <a onClick="popup('http://www.biomedcentral.com/1471-2105/14/99/mathml/M67','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/14/99/mathml/M67">View MathML</a>. The second step of the algorithm consists of estimating the unknown model parameters (M, A and σ2) given<a onClick="popup('http://www.biomedcentral.com/1471-2105/14/99/mathml/M68','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/14/99/mathml/M68">View MathML</a>. The estimated posterior distributions of the unknown model parameters are given in Additional file 1: Figure S5 and validate the proposed Bayesian model.

The burn-in period and number of Gibbs samples were determined using quantitative methods described in the Additional file 1: Section “Convergence diagnosis”.

Comparison to other methods

The performance of the proposed uBLU algorithm is compared with other existing factor decomposition methods including PCA, NMF, BFRM and GB-GMF by using the following criteria, which are common measures used to compare factor analysis algorithms,

• the factor mean square errors (MSE)

<a onClick="popup('http://www.biomedcentral.com/1471-2105/14/99/mathml/M69','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/14/99/mathml/M69">View MathML</a>

where<a onClick="popup('http://www.biomedcentral.com/1471-2105/14/99/mathml/M70','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/14/99/mathml/M70">View MathML</a> is the estimated r-th factor loading vector,

• the global MSE of factor scores

<a onClick="popup('http://www.biomedcentral.com/1471-2105/14/99/mathml/M71','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/14/99/mathml/M71">View MathML</a>

where <a onClick="popup('http://www.biomedcentral.com/1471-2105/14/99/mathml/M72','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/14/99/mathml/M72">View MathML</a> is the estimated proportion of the r-th factor in the i-th sample,

• the reconstruction error (RE)

<a onClick="popup('http://www.biomedcentral.com/1471-2105/14/99/mathml/M73','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/14/99/mathml/M73">View MathML</a>

(21)

• where<a onClick="popup('http://www.biomedcentral.com/1471-2105/14/99/mathml/M74','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/14/99/mathml/M74">View MathML</a> is the estimate of yi,

• the spectral angle distance (SAD) between mr and its estimate<a onClick="popup('http://www.biomedcentral.com/1471-2105/14/99/mathml/M75','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/14/99/mathml/M75">View MathML</a> for each factor r = 1, …, R

<a onClick="popup('http://www.biomedcentral.com/1471-2105/14/99/mathml/M76','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/14/99/mathml/M76">View MathML</a>

where arccos(·) is the inverse cosine function,

• the global spectral angle distance (GSAD) between yi (the i-th observation vector) and <a onClick="popup('http://www.biomedcentral.com/1471-2105/14/99/mathml/M77','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/14/99/mathml/M77">View MathML</a> (its estimate)

<a onClick="popup('http://www.biomedcentral.com/1471-2105/14/99/mathml/M78','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/14/99/mathml/M78">View MathML</a>

• the computational time.

The proposed uBLU algorithm, the PCA, NMF and GB-GMF methods were implemented in Matlab 7.8.0 (R2009a). The BFRM software (version 2.0) was downloaded from [22] and implemented with default values for the parameters. All methods were implemented on an Intel(R) Core(TM)2 Duo processor.

Simulation results are reported in Tables 2, 3, 4 and 5. Note that the positivity and sum-to-one constraints that are enforced on the data for the proposed uBLU algorithm avoid the scale ambiguity inherent to any factor decomposition problem. Conversely, for the other factor decomposition methods (PCA, NMF, BFRM and GB-GMF), if {M, A} is an admissible solution, {MB, BTA} is also admissible for any scaling and permutation matrix B. Hence a re-scaling is required to identify appropriate permutations before computing MSEs and GMSEs. Moreover, when PCA, NMF, BFRM and GB-GMF methods are run for R = 4, we only considered the 3 factors yielding the 3 smallest SADs values.

Table 2. Simulation results for dataset<a onClick="popup('http://www.biomedcentral.com/1471-2105/14/99/mathml/M79','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/14/99/mathml/M79">View MathML</a>

Table 3. Simulation results for dataset<a onClick="popup('http://www.biomedcentral.com/1471-2105/14/99/mathml/M89','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/14/99/mathml/M89">View MathML</a>

Table 4. Simulation results for dataset<a onClick="popup('http://www.biomedcentral.com/1471-2105/14/99/mathml/M99','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/14/99/mathml/M99">View MathML</a>

Table 5. Simulation results for dataset<a onClick="popup('http://www.biomedcentral.com/1471-2105/14/99/mathml/M109','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/14/99/mathml/M109">View MathML</a>

These results show that the uBLU method is more flexible since it provides better unmixing performance across all of the considered synthetic datasets <a onClick="popup('http://www.biomedcentral.com/1471-2105/14/99/mathml/M119','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/14/99/mathml/M119">View MathML</a> as compared to other existing factorization methods (PCA, NMF, BFRM and GB-GMF). Moreover, uBLU has the following advantages: i) it is fully unsupervised and does not require the number of factors to be specified as a prior knowledge, ii) due to the constraints, the factors and factor scores are estimated without scale ambiguity. The disadvantage is the execution time: uBLU requires more computation due to the Gibbs sampling.

Evaluation on gene expression data

Here the proposed algorithm is illustrated on a real time-evolving gene expression data from recent viral challenge studies on influenza A/H3N2/Wisconsin. The data are available at GEO, accession number GSE30550.

Details on data collection

We briefly describe the dataset. For more details the reader is referred to [14,18]. H3N2 dataset consists of the gene expression levels of N = 267 Affymetrix chips collected on 17 healthy human volunteers experimentally infected with influenza A/Wisconsin/67/2005 (H3N2). A clinical symptom score was assigned to each sample by clinicians who participated in the study. Nine of the 17 subjects (those labeled Z01, Z05, Z06, Z07, Z08, Z10, Z12, Z13, and Z15 in Figure 1c) became clinically ill during the study. These labels are only used as ground truth to quantify performance and are not available to the uBLU algorithm. The challenge consists of inoculating intranasally a dose of 106 TCID50 Influenza A manufactured and processed under current good manufacturing practices (cGMP) by Baxter BioScience. Peripheral blood microarray analysis was performed at multiple time instants corresponding to baseline (24 hours prior to inoculation with virus), then at 8 hour intervals for the initial 120 hours and then 24 hours for two further days. Each sample consisted of over G = 12000 gene expression values after standard microarray data normalization with RMA using the custom brain array cdf [14]. No other preprocessing was applied prior to running the five unsupervised methods (uBLU, PCA, NMF, BFRM, and GB-GMF).

thumbnailFigure 1. Experimental results on the H3N2 viral challenge dataset of gene expression profiles. (a) Estimated posterior distribution of the number of factors R. (b) Factor loadings ranked by decreasing dominance. (c) Heatmap of the factor scores of the inflammatory component clearly separates symptomatic subjects (bottom 9 rows) and the time course of their molecular inflammatory response. The five black colored pixels indicate samples that were not assayed.

Application of the proposed uBLU algorithm

The uBLU algorithm was run with Nmc = 10000 Monte Carlo iterations, including a burn-in period of Nbi = 1000 iterations. uBLU allows the posterior distribution of the number of factors R, depicted in Figure 1a, to be estimated. The results show that the MAP estimate of the number of factors is<a onClick="popup('http://www.biomedcentral.com/1471-2105/14/99/mathml/M120','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/14/99/mathml/M120">View MathML</a> (more than 90% of the generated Gibbs samples of the number of factors were equal to 4).

Figure 2 shows the reconstruction error RE(t) as a function of the number of iterations (t = 1, …). The reconstruction errors are computed from the observed gene expression data matrix and the estimates of the factor and factor score matrix M and A at each iteration. Figure 2 also indicates that the number of burn-in and Monte Carlo samples Nbi = 1000 and Nmc = 10000 are sufficient.

thumbnailFigure 2. Reconstruction error and estimated number of factors as a function of the number of iterations (H3N2 challenge data). Top: Reconstruction error (RE(t)) computed from the observation matrix Y and the estimated matrices M(t) and A(t) as a function of the iteration index t. Bottom: Estimated number of factors R(t) as a function of the iteration number t.

The different factors are depicted in Figure 1b where the G genes have been reordered so that the dominant genes are grouped together in each factor. Factors are then ordered with respect to their maximum loading. Specifically, the k-th sharp peak in the figure occurs at the gene index that has maximal loading in factor k. Genes to the right of this dominant gene up to the (k + 1)-st peak also dominate in this k-th factor, but at a lower degree. uBLU identifies a strong factor (the first factor, in red) by virtue of its significantly larger proportion of highly dominant genes. Many of the genes in this strong factor are recognizable as immune response genes that regulate pattern recognition, interferon, and inflammation pathways in respiratory viral response. A very similar factor was found in a different analysis [14,18] of this dataset and here we call it the “inflammatory component”.

The factor scores corresponding to this inflammatory component are shown in Figure 1c, where they are rendered as an image whose columns (respectively rows) index the subjects (respectively the different time sampling instants). Figure 1c shows that uBLU clearly separates the samples of subjects exhibiting symptoms (associated with the last 9 rows) from those who remain asymptomatic (associated with the first 8 rows), when the estimated number of factors is<a onClick="popup('http://www.biomedcentral.com/1471-2105/14/99/mathml/M121','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/14/99/mathml/M121">View MathML</a>. Moreover, this symptom factor can be used to segment the data matrix into 3 states: pre-onset-symptomatic (before significant symptoms occur), post-onset-symptomatic and asymptomatic.

Furthermore, this inflammatory factor identified by the proposed uBLU algorithm is most highly represented in those samples associated with acute flu symptoms, as measured by modified Jackson scores (see [14], Figure 1B). The dominant gene contributors to this inflammatory component correspond to well-known transcription factors controlling immune response, inflammatory response and antigen presentation – see Table 6. The reader is referred to [14,18] for more details on clinical determination of symptom scores and biological significance of the inflammatory component genes.

Table 6. NCI-curated pathway associations of group of genes contributing to uBLU inflammatory component

For comparison we applied a supervised version of the proposed uBLU algorithm to the H3N2 dataset. This was implemented by setting the number of factors to R = 4 and using the algorithm [13] to jointly estimate M and A. The inflammatory component found by the supervised algorithm was virtually identical to the one found by the proposed algorithm (uBLU) that automatically selects R = 4.

Comparison to other methods

The uBLU algorithm is compared with four matrix factorization algorithms, i.e. PCA, NMF, BFRM and GB-GMF methods.

Figure 3 depicts the different factors, ordered so that the inflammatory group is the leftmost one (in red). The factor loadings obtained with NMF or PCA reveal the inflammatory component. However, there are fewer highly dominant genes in the NMF and PCA loadings for this factor as compared to uBLU. The BFRM and GB-GMF methods found four pathways, several overlapping with those of uBLU, NMF and PCA.

thumbnailFigure 3. Factor loadings ranked by decreasing dominance for H3N2 challenge data. uBLU shows a particularly strong component (Figure 1b), the group 1, that corresponds to the well-known inflammatory pathway. NMF and PCA algorithms also reveal an inflammatory component, but it includes fewer relevant genes than uBLU. See Figure 4 for the corresponding factor scores.

thumbnailFigure 4. Heatmaps of the factor scores of the inflammatory component for H3N2 challenge data. The inflammatory factor determined by the proposed uBLU method (a) shows higher contrast between symptomatic and asymptomatic subjects than the other methods. The five black colored pixels of the heatmaps indicate samples that were not assayed.

The factor scores of the five matrix factorization methods corresponding to the inflammatory component are depicted in Figure 4. This figure shows that the uBLU and the NMF methods are better able to attain a high contrast separation between the acutely symptomatic samples and the other samples. This is confirmed by the evaluation of the Fisher criteria (22) between these two regions (see Table 7). Indeed, denote by (μpos, σ2pos) the empirical mean and variance of the scores associated with the Npos samples in the acute symptomatic state (bright colored samples in the lower right rectangle of Figure 1c). Denote by <a onClick="popup('http://www.biomedcentral.com/1471-2105/14/99/mathml/M122','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/14/99/mathml/M122">View MathML</a> the same parameters for the remaining samples. The Fisher linear discriminant measure ([23], p. 119) is defined as

<a onClick="popup('http://www.biomedcentral.com/1471-2105/14/99/mathml/M123','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/14/99/mathml/M123">View MathML</a>

(22)

Table 7. Simulation results for real H3N2 dataset

To compare the biological relevance of the inflammatory genes found by uBLU to those found by the other methods we performed gene enrichment analysis (GEA). Here we only report GEA comparisons between uBLU and NMF. Tables 6 and 8 show the pathway enrichment associated with the top 200 genes found by uBLU and NMF, respectively, using NCI pathway interaction database (http://pid.nci.nih.gov webcite). The uBLU genes are significantly better associated with the NCI-curated pathways than the NMF genes. In particular, the two most enriched pathways, IFN-gamma and PDGFR beta signaling, associated with the uBLU genes have much higher statistical significance (lower p-value) than any of the pathways associated with NMF.

Table 8. NCI-curated pathway associations of group of genes contributing to NMF inflammatory component

Figure 5 shows how the factor scores of the dominant factor can be used as features to cluster samples. Euclidean multidimensional scaling (MDS) [24] is used to map the factor score vector for each sample as a coordinate in the plane. Each sample is embedded with a color and a size, denoting the state of the subject (asymptomatic subjects in blue, symptomatic subjects in red) and the time stamp, respectively. These figures show that uBLU can separate sick and healthy subjects, as well as or better than NMF, BFRM and GB-GMF.

thumbnailFigure 5. Chip clouds after demixing for H3N2 challenge data. These figures show the scatter of the four dimensional factor score vectors (projected onto the plane using MDS) for each algorithm that was compared to uBLU. uBLU, NMF and BFRM obtain a clean separation of samples of symptomatic (red points) and asymptomatic (blue points) subjects whereas the separation is less clear with PCA. In these scatter plots the size of a point is proportional to the time at which the sample was taken during challenge study.

One can conclude from these comparisons that, when applied on the H3N2 dataset, the proposed uBLU algorithm outperforms PCA, NMF, BFRM, and GB-GMF algorithms in terms of finding genes with higher pathway enrichment and achieving higher contrast of the acute symptom states.

The computational times required by the five considered matrix factorization methods, including the proposed uBLU algorithm, when applied to this real dataset, are reported in Table 7. The GB-GMF algorithm is slightly faster than the proposed algorithm but does not identify the inflammatory component or achieve good contrast of the acute symptom states in the H3N2 challenge study.

Conclusions

This paper proposes a new Bayesian unmixing algorithm for discovering signatures in high dimensional biological data, and specifically for gene expression microarrays. An interesting property of the proposed algorithm is that it provides positive factor loadings to ensure positivity as well as sum-to-one constraints for the factor scores. The advantages of these constraints are that they lead to better discrimination between sick and healthy individuals, and they recover the inflammatory genes in a unique factor, the inflammatory component. The proposed algorithm is fully unsupervised in the sense that it does not depend on any labeling of the samples and that it can infer the number of factors directly from the observation data matrix. Finally, as any Bayesian algorithm, the Monte Carlo-based procedure investigated in this study provides point estimates as well as confidence intervals for the unknown parameters, contrary to many existing factor decomposition methods such as PCA or NMF.

Simulation results performed on synthetic and real data demonstrated significant improvements. Indeed, when applied to real time-evolving gene expression datasets, the uBLU algorithm revealed an inflammatory factor with higher contrast between subjects who would become symptomatic from those who would remain asymptomatic (as determined by comparing to ground truth clinical labels).

In this study, the time samples were modeled as independent. Future works include extensions of the proposed model to account for time dependency between samples.

Appendix A: Gibbs sampler

This appendix provides more details about the Gibbs sampler strategy to generate samples<a onClick="popup('http://www.biomedcentral.com/1471-2105/14/99/mathml/M124','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/14/99/mathml/M124">View MathML</a> distributed according to the joint distribution<a onClick="popup('http://www.biomedcentral.com/1471-2105/14/99/mathml/M125','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/14/99/mathml/M125">View MathML</a> (the reader is referred to [25] for more details about the Gibbs sampler and MCMC methods). This joint distribution can be obtained by integrating out the hyperparameter γ from f(Θ, γ|Y) defined in (18) and can be written

<a onClick="popup('http://www.biomedcentral.com/1471-2105/14/99/mathml/M126','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/14/99/mathml/M126">View MathML</a>

(23)

where the dimensions of the matrices M, T, and A depend on the unknown number of factors R and the priors have been defined in the Section “Methods”.

The different steps of the Gibbs sampler are detailed below.

Inference of the number of factors

The proposed unsupervised algorithm includes a birth/death process for inferring the number of factors R, i.e., it generates samples R in addition to M and A. More precisely, at iteration t of the algorithm, a birth, death or switch move is randomly chosen with probabilities <a onClick="popup('http://www.biomedcentral.com/1471-2105/14/99/mathml/M127','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/14/99/mathml/M127">View MathML</a>, <a onClick="popup('http://www.biomedcentral.com/1471-2105/14/99/mathml/M128','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/14/99/mathml/M128">View MathML</a> and <a onClick="popup('http://www.biomedcentral.com/1471-2105/14/99/mathml/M129','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/14/99/mathml/M129">View MathML</a>. The birth and death moves consist of increasing or decreasing by 1 the number R of factors using a reversible jump step (see [21] for more details), whereas the switch move does not change the dimension of R and requires the use of a Metropolis-Hastings acceptance procedure. Let consider a move, at iteration index t, from the state {M(t), A(t), R(t)} to the new state {M, A, R}. The birth, death and switch moves are defined as follows, similar to those used in [26] (Algorithms 3, 4 and 5).

Birth move: When a birth move is proposed, a new signature m is randomly generated to build M = [M(t), m]. The new corresponding space is checked so that the signatures are sufficiently distinct and separate from one another. Then, a new factor score coefficient is drawn, for each vector ai (i = 1, …, N), from a Beta distribution <a onClick="popup('http://www.biomedcentral.com/1471-2105/14/99/mathml/M130','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/14/99/mathml/M130">View MathML</a>, and the new factor score matrix, denoted as A, is re-scaled to sum to one.

Death move: When a death move is proposed, one of the factors of M(t), and its corresponding factor score coefficients, are randomly removed. The remaining factor scores are re-scaled to ensure the sum-to-one constraint.

Switch move: When a switch move is proposed, a signature m is randomly chosen and replaced with another signature randomly generated. If the new signature is too close to another, its corresponding factor scores are proportionately distributed among its closest factors. Indeed, the switch move consists of creating a new signature (birth move) and deleting another one (death move) in a faster single step.

Each move is then accepted or rejected according to an empirical acceptance probability: the likelihood ratio between the actual state and the proposed new state. The factor matrix M, the factor score matrix A and the noise variance σ2 are then updated, conditionally upon the number of factors R, using the following Gibbs steps.

Generation of samples according to f(T|A, σ2, R, Y)

Sampling from the joint conditional f(T|A, σ2, R,Y) is achieved by updating each column of T using Gibbs moves. Let denote Tr the matrix T whose r-th column has been removed. The posterior distribution of tr is the following truncated multivariate Gaussian distribution (MGD)

<a onClick="popup('http://www.biomedcentral.com/1471-2105/14/99/mathml/M131','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/14/99/mathml/M131">View MathML</a>

(24)

where

<a onClick="popup('http://www.biomedcentral.com/1471-2105/14/99/mathml/M132','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/14/99/mathml/M132">View MathML</a>

(25)

For more details on how we generate realizations from this truncated distribution, see [13].

Generation of samples according to f(a1:R − 1, i|T, σ2, R, Y)

Straightforward computations lead to the posterior distribution of each element of a1:R − 1, i

<a onClick="popup('http://www.biomedcentral.com/1471-2105/14/99/mathml/M133','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/14/99/mathml/M133">View MathML</a>

(26)

where

<a onClick="popup('http://www.biomedcentral.com/1471-2105/14/99/mathml/M134','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/14/99/mathml/M134">View MathML</a>

(27)

<a onClick="popup('http://www.biomedcentral.com/1471-2105/14/99/mathml/M135','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/14/99/mathml/M135">View MathML</a> and MR denotes the factor loading matrix M whose R-th column has been removed. Equation (26) shows that the factor score distribution is an MGD truncated on the simplex <a onClick="popup('http://www.biomedcentral.com/1471-2105/14/99/mathml/M136','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/14/99/mathml/M136">View MathML</a> defined in (12).

Generation of samples according to f(σ2|T, A, R, Y)

Using (14) and (16), one can show that the conditional distribution f(σ2|M, A, Y) is the following inverse-Gamma distribution

<a onClick="popup('http://www.biomedcentral.com/1471-2105/14/99/mathml/M137','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2105/14/99/mathml/M137">View MathML</a>

(28)

Appendix B: Contribution of each of uBLU’s constraints

To illustrate the advantage of enforcing non-negativity and sum-to-one constraints on the factors and on the factor scores, as detailed in the Methods section, we evaluated the effect of successively stripping out these constraints from uBLU. In particular we implemented uBLU under the following conditions: i) without any constraints, ii) with only the positivity constraints on the factors and the scores, iii) with only the sum-to-one constraint on the scores, and iv) with both positivity and sum-to-one constraint on factors and scores as proposed in (5).

Figures 6 display heatmaps of the factor scores of the inflammatory component. The segmentation into two main regions (post-symptomatic samples and asymptomatic samples) becomes apparent only when the sum-to-one constraint is enforced on the scores. To quantify the benefit that is visible in Figure 6 we performed a GEA analysis, reported in Table 9, on the top 200 genes found in each of the inflammatory components found by uBLU implemented with no constraints, positivity constraints, sum-to-one constraints, and both constraints. The table shows that both constraints are necessary to obtain the best enrichment scores (lowest possible p-values).

thumbnailFigure 6. Contribution of each constraint on the scores of the inflammatory factor (H3N2 challenge data). The five black colored pixels of the heatmaps indicate samples that were not assayed. Note that when only the sum-to-one constraint is applied, non-negativity is not guaranteed. However, for this dataset the sum-to-one factor scores turn out to take on non-negative values for the inflammatory factor (but not for the other factors).

Table 9. Contribution of each of uBLU’s constraints

Competing interests

The authors declare that they have no competing interests.

Authors’ contributions

CB, ND, JYT and AH performed the statistical analysis. GG and AZ designed the Flu challenge experiment that generated the data used to compare the methods. All authors contributed to the manuscript and approved the final version.

Acknowledgements

This work was supported in part by DARPA under the PHD program (DARPA-N66001-09-C-2082). The views, opinions, and findings contained in this article are those of the authors and should not be interpreted as representing the official views or policies, either expressed or implied, of the Defense Advanced Research Projects Agency or the Department of Defense. Approved for Public Release, Distribution Unlimited.

References

  1. Carvalho CM, Chang J, Lucas JE, Nevins JR, Wang Q, West M: High-dimensional sparse factor modelling: applications in gene expression genomics.

    J Am Stat Assoc 2008, 103(484):1438-1456. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  2. Paisley J, Carin L: Nonparametric factor analysis with beta process priors. In Proc 26th Annual Int Conf on Machine Learning. ICML 2009. Montreal, Quebec, Canada; 2009:777-784. OpenURL

  3. Chen B, Chen M, Paisley J, Zaas A, Woods C, Ginsburg GS, Hero AO, Lucas J, Dunson D, Carin L: Bayesian inference of the number of factors in gene-expression analysis: application to human virus challenge studies.

    BMC Bioinformatics 2010, 11:552. PubMed Abstract | BioMed Central Full Text | PubMed Central Full Text OpenURL

  4. West M: Bayesian Factor regression models in the “Large p, Small n” paradigm. In Bayesian Statistics 7. Oxford University Press; 2003:723-732. OpenURL

  5. Yeung KY, Ruzzo WL: Principal component analysis for clustering gene expression data.

    Bioinformatics 2001, 17(9):763-774. PubMed Abstract | Publisher Full Text OpenURL

  6. Nascimento JM, Bioucas-Dias JM: Vertex component analysis: A fast algorithm to unmix hyperspectral data.

    IEEE Trans Geosci Remote Sensing 2005, 43(4):898-910. OpenURL

  7. Lee DD, Seung HS: Algorithms for non-negative matrix factorization.

    Proc Neural Info Process Syst 2000, 13:556-562. OpenURL

  8. Fogel P, Young SS, Hawkins DM, Ledirac N: Inferential, robust non-negative matrix factorization analysis of microarray data.

    Bioinformatics 2007, 23:44-49. PubMed Abstract | Publisher Full Text OpenURL

  9. McLachlan GJ, Bean RW, Peel D: A mixture model-based approach to the clustering of microarray expression data.

    Bioinformatics 2002, 18(3):413-422. PubMed Abstract | Publisher Full Text OpenURL

  10. Baek J, McLachlan GJ: Mixtures of common t-factor analyzers for clustering high-dimensional microarray data.

    Bioinformatics 2011, 27:1269-1276. PubMed Abstract | Publisher Full Text OpenURL

  11. Moloshok TD, Klevecz RR, Grant JD, Manion FJ, Speier WF, Ochs MF: Application of Bayesian decomposition for analysing microarray data.

    Bioinformatics 2002, 18:566-575. PubMed Abstract | Publisher Full Text OpenURL

  12. Fertig EJ, Ding J, Favorov AV, Parmigiani G, Ochs MF: CoGAPS: an R/C++ package to identify patterns and biological process activity in transcriptomic data.

    Bioinformatics 2010, 26:2792-2793. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  13. Dobigeon N, Moussaoui S, Coulon M, Tourneret JY, Hero AO: Joint Bayesian endmember extraction and linear unmixing for hyperspectral imagery.

    IEEE Trans Signal Process 2009, 57(11):4355-4368. OpenURL

  14. Huang Y, Zaas AK, Rao A, Dobigeon N, Woolf PJ, Veldman T, Oien NC, McClain MT, Varkey JB, Nicholson B, Carin L, Kingsmore S, Woods CW, Ginsburg GS, Hero A: Temporal dynamics of host molecular responses differentiate symptomatic and asymptomatic influenza A infection.

    PLoS Genet 2011, 8(7):e1002234. OpenURL

  15. Hyvärinen A, Karhunen J, Oja E: Independent Component Analysis. New York: John Wiley; 2001. OpenURL

  16. Dueck D, Morris QD, Frey BJ: Multi-way clustering of microarray data using probabilistic sparse matrix factorization.

    Bioinformatics 2005, 21:144-151. PubMed Abstract | Publisher Full Text OpenURL

  17. Nikulin V, Huang TH, Ng SK, Rathnayake S, McLachlan GJ: A very fast algorithm for matrix factorization.

    Stat Probability Lett 2011, 81(7):773-782. Publisher Full Text OpenURL

  18. Zaas AK, Chen M, Varkey J, Veldman T, Hero AO, Lucas J, Huang Y, Turner R, Gilbert A, Lambkin-Williams R, Øien NC, Nicholson B, Kingsmore S, Carin L, Woods CW, Ginsburg GS: Gene expression signatures diagnose influenza and other symptomatic respiratory viral infections in humans.

    Cell Host Microbe 2009, 6(3):207-217.

    [http://www.ncbi.nlm.nih.gov/pubmed/19664979 webcite]

    PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  19. Winter ME: N-FINDR: an algorithm for fast autonomous spectral end-member determination in hyperspectral data.

    Imaging Spectrometry V Proc SPIE 3753 1999, 266-275. OpenURL

  20. Gilks WR, Richardson S, Spiegelhalter DJ: Markov Chain Monte Carlo in Practice. London: Chapman and Hall; 1996.

    (ISBN: 0-412-05551-1)

    OpenURL

  21. Green PJ: Reversible jump Markov chain Monte Carlo computation and Bayesian model determination.

    Biometrika 1995, 82(4):711-732. Publisher Full Text OpenURL

  22. BFRM Software: bayesian factor regression modelling [http://www.stat.duke.edu/research/software/west/bfrm/download.html webcite]

  23. Duda RO, Hart PE, Stork DG: Pattern Classification. New York: Wiley-Interscience; 2000. OpenURL

  24. Cox TF, Cox MAA: Multidimensional Scaling. London: Chapman and Hall; 1994. OpenURL

  25. Robert CP, Casella G: Monte Carlo Statistical Methods. New York: Springer-Verlag; 1999. OpenURL

  26. Dobigeon N, Tourneret JY, Chang CI: Semi-supervised linear spectral unmixing using a hierarchical Bayesian model for hyperspectral imagery.

    IEEE Trans Signal Process 2008, 56(7):2684-2695. OpenURL