Generative models on discrete state-spaces have a wide range of potential applications, particularly in the domain of natural sciences. In continuous state-spaces, controllable and flexible generation of samples with desired properties has been realized using guidance on diffusion and flow models. However, these guidance approaches are not readily amenable to discrete state-space models. Consequently, we introduce a general and principled method for applying guidance on such models. Our method depends on leveraging continuous-time Markov processes on discrete state-spaces, which unlocks computational tractability for sampling from a desired guided distribution. We demonstrate the utility of our approach, Discrete Guidance, on a range of applications including guided generation of small-molecules, DNA sequences and protein sequences.

SUMMARY: The design of two overlapping genes in a microbial genome is an emerging technique for adding more reliable control mechanisms in engineered organisms for increased stability. The design of functional overlapping gene pairs is a challenging procedure, and computational design tools are used to improve the efficiency to deploy successful designs in genetically engineered systems. GENTANGLE (Gene Tuples ArraNGed in overLapping Elements) is a high-performance containerized pipeline for the computational design of two overlapping genes translated in different reading frames of the genome. This new software package can be used to design and test gene entanglements for microbial engineering projects using arbitrary sets of user-specified gene pairs. AVAILABILITY AND IMPLEMENTATION: The GENTANGLE source code and its submodules are freely available on GitHub at https://github.com/BiosecSFA/gentangle. The DATANGLE (DATA for genTANGLE) repository contains related data and results and is freely available on GitHub at https://github.com/BiosecSFA/datangle. The GENTANGLE container is freely available on Singularity Cloud Library at https://cloud.sylabs.io/library/khyox/gentangle/gentangle.sif. The GENTANGLE repository wiki (https://github.com/BiosecSFA/gentangle/wiki), website (https://biosecsfa.github.io/gentangle/), and user manual contain detailed instructions on how to use the different components of software and data, including examples and reproducing the results. The code is licensed under the GNU Affero General Public License version 3 (https://www.gnu.org/licenses/agpl.html).

Data-driven design is making headway into a number of application areas, including protein, small-molecule, and materials engineering. The design goal is to construct an object with desired properties, such as a protein that binds to a therapeutic target, or a superconducting material with a higher critical temperature than previously observed. To that end, costly experimental measurements are being replaced with calls to high-capacity regression models trained on labeled data, which can be leveraged in an in silico search for design candidates. However, the design goal necessitates moving into regions of the design space beyond where such models were trained. Therefore, one can ask: should the regression model be altered as the design algorithm explores the design space, in the absence of new data? Herein, we answer this question in the affirmative. In particular, we (i) formalize the data-driven design problem as a non-zero-sum game, (ii) develop a principled strategy for retraining the regression model as the design algorithm proceeds---what we refer to as autofocusing, and (iii) demonstrate the promise of autofocusing empirically.

We present a method for design problems wherein the goal is to maximize or specify the value of one or more properties of interest (e.g. maximizing the fluorescence of a protein). We assume access to black box, stochastic “oracle" predictive functions, each of which maps from design space to a distribution over properties of interest. Because many state-of-the-art predictive models are known to suffer from pathologies, especially for data far from the training distribution, the problem becomes different from directly optimizing the oracles. Herein, we propose a method to solve this problem that uses model-based adaptive sampling to estimate a distribution over the design space, conditioned on the desired properties.

Variational autoencoders (VAE) are a powerful and widely-used class of models to learn complex data distributions in an unsupervised fashion. One important limitation of VAEs is the prior assumption that latent sample representations are independent and identically distributed. However, for many important datasets, such as time-series of images, this assumption is too strong: accounting for covariances between samples, such as those in time, can yield to a more appropriate model specification and improve performance in downstream tasks. In this work, we introduce a new model, the Gaussian Process (GP) Prior Variational Autoencoder (GPPVAE), to specifically address this issue. The GPPVAE aims to combine the power of VAEs with the ability to model correlations afforded by GP priors. To achieve efficient inference in this new class of models, we leverage structure in the covariance matrix, and introduce a new stochastic backpropagation strategy that allows for computing stochastic gradients in a distributed and low-memory fashion. We show that our method outperforms conditional VAEs (CVAEs) and an adaptation of standard VAEs in two image data applications.

MOTIVATION: Set-based variance component tests have been identified as a way to increase power in association studies by aggregating weak individual effects. However, the choice of test statistic has been largely ignored even though it may play an important role in obtaining optimal power. We compared a standard statistical test-a score test-with a recently developed likelihood ratio (LR) test. Further, when correction for hidden structure is needed, or gene-gene interactions are sought, state-of-the art algorithms for both the score and LR tests can be computationally impractical. Thus we develop new computationally efficient methods. RESULTS: After reviewing theoretical differences in performance between the score and LR tests, we find empirically on real data that the LR test generally has more power. In particular, on 15 of 17 real datasets, the LR test yielded at least as many associations as the score test-up to 23 more associations-whereas the score test yielded at most one more association than the LR test in the two remaining datasets. On synthetic data, we find that the LR test yielded up to 12% more associations, consistent with our results on real data, but also observe a regime of extremely small signal where the score test yielded up to 25% more associations than the LR test, consistent with theory. Finally, our computational speedups now enable (i) efficient LR testing when the background kernel is full rank, and (ii) efficient score testing when the background kernel changes with each test, as for gene-gene interaction tests. The latter yielded a factor of 2000 speedup on a cohort of size 13 500. AVAILABILITY: Software available at http://research.microsoft.com/en-us/um/redmond/projects/MSCompBio/Fastlmm/. CONTACT: [email protected] SUPPLEMENTARY INFORMATION: Supplementary data are available at Bioinformatics online.

MOTIVATION: Approaches for testing sets of variants, such as a set of rare or common variants within a gene or pathway, for association with complex traits are important. In particular, set tests allow for aggregation of weak signal within a set, can capture interplay among variants and reduce the burden of multiple hypothesis testing. Until now, these approaches did not address confounding by family relatedness and population structure, a problem that is becoming more important as larger datasets are used to increase power. RESULTS: We introduce a new approach for set tests that handles confounders. Our model is based on the linear mixed model and uses two random effects-one to capture the set association signal and one to capture confounders. We also introduce a computational speedup for two random-effects models that makes this approach feasible even for extremely large cohorts. Using this model with both the likelihood ratio test and score test, we find that the former yields more power while controlling type I error. Application of our approach to richly structured Genetic Analysis Workshop 14 data demonstrates that our method successfully corrects for population structure and family relatedness, whereas application of our method to a 15 000 individual Crohn's disease case-control cohort demonstrates that it additionally recovers genes not recoverable by univariate analysis. AVAILABILITY: A Python-based library implementing our approach is available at http://mscompbio.codeplex.com.

High-resolution HLA typing plays a central role in many areas of immunology, such as in identifying immunogenetic risk factors for disease, in studying how the genomes of pathogens evolve in response to immune selection pressures, and also in vaccine design, where identification of HLA-restricted epitopes may be used to guide the selection of vaccine immunogens. Perhaps one of the most immediate applications is in direct medical decisions concerning the matching of stem cell transplant donors to unrelated recipients. However, high-resolution HLA typing is frequently unavailable due to its high cost or the inability to re-type historical data. In this paper, we introduce and evaluate a method for statistical, in silico refinement of ambiguous and/or low-resolution HLA data. Our method, which requires an independent, high-resolution training data set drawn from the same population as the data to be refined, uses linkage disequilibrium in HLA haplotypes as well as four-digit allele frequency data to probabilistically refine HLA typings. Central to our approach is the use of haplotype inference. We introduce new methodology to this area, improving upon the Expectation-Maximization (EM)-based approaches currently used within the HLA community. Our improvements are achieved by using a parsimonious parameterization for haplotype distributions and by smoothing the maximum likelihood (ML) solution. These improvements make it possible to scale the refinement to a larger number of alleles and loci in a more computationally efficient and stable manner. We also show how to augment our method in order to incorporate ethnicity information (as HLA allele distributions vary widely according to race/ethnicity as well as geographic area), and demonstrate the potential utility of this experimentally. A tool based on our approach is freely available for research purposes at http://microsoft.com/science.

MOTIVATION: There is a pressing need for improved proteomic screening methods allowing for earlier diagnosis of disease, systematic monitoring of physiological responses and the uncovering of fundamental mechanisms of drug action. The combined platform of LC-MS (Liquid-Chromatography-Mass-Spectrometry) has shown promise in moving toward a solution in these areas. In this paper we present a technique for discovering differences in protein signal between two classes of samples of LC-MS serum proteomic data without use of tandem mass spectrometry, gels or labeling. This method works on data from a lower-precision MS instrument, the type routinely used by and available to the community at large today. We test our technique on a controlled (spike-in) but realistic (serum biomarker discovery) experiment which is therefore verifiable. We also develop a new method for helping to assess the difficulty of a given spike-in problem. Lastly, we show that the problem of class prediction, sometimes mistaken as a solution to biomarker discovery, is actually a much simpler problem. RESULTS: Using precision-recall curves with experimentally extracted ground truth, we show that (1) our technique has good performance using seven replicates from each class, (2) performance degrades with decreasing number of replicates, (3) the signal that we are teasing out is not trivially available (i.e. the differences are not so large that the task is easy). Lastly, we easily obtain perfect classification results for data in which the problem of extracting differences does not produce absolutely perfect results. This emphasizes the different nature of the two problems and also their relative difficulties. AVAILABILITY: Our data are publicly available as a benchmark for further studies of this nature at http://www.cs.toronto.edu/~jenn/LCMS

The identification of T cell epitopes and their HLA (human leukocyte antigen) restrictions is important for applications such as the design of cellular vaccines for HIV. Traditional methods for such identification are costly and time-consuming. Recently, a more expeditious laboratory technique using ELISpot assays has been developed that allows for rapid screening of specific responses. However, this assay does not directly provide information concerning the HLA restriction of a response, a critical piece of information for vaccine design. Thus, we introduce, apply, and validate a statistical model for identifying HLA-restricted epitopes from ELISpot data. By looking at patterns across a broad range of donors, in conjunction with our statistical model, we can determine (probabilistically) which of the HLA alleles are likely to be responsible for the observed reactivities. Additionally, we can provide a good estimate of the number of false positives generated by our analysis (i.e., the false discovery rate). This model allows us to learn about new HLA-restricted epitopes from ELISpot data in an efficient, cost-effective, and high-throughput manner. We applied our approach to data from donors infected with HIV and identified many potential new HLA restrictions. Among 134 such predictions, six were confirmed in the lab and the remainder could not be ruled as invalid. These results shed light on the extent of HLA class I promiscuity, which has significant implications for the understanding of HLA class I antigen presentation and vaccine development.

We present a hierarchical Bayesian model for sets of related, but different, classes of time series data. Our model performs alignment simultaneously across all classes, while detecting and characterizing class-specific differences. During inference the model produces, for each class, a distribution over a canonical representation of the class. These class-specific canonical representations are automatically aligned to one another -- preserving common sub-structures, and highlighting differences. We apply our model to compare and contrast solenoid valve current data, and also, liquid-chromatography-ultraviolet-diode array data from a study of the plant Arabidopsis thaliana. 1 Aligning Time Series From Different Classes Many practical problems over a wide range of domains require synthesizing information from several noisy examples of one or more categories in order to build a model which captures common structure and also learns the patterns of variability between categories. In time series analysis, these modeling goals manifest themselves in the tasks of alignment and difference detection. These tasks have diverse applicability, spanning speech & music processing, equipment & industrial plant diagnosis/monitoring, and analysis of biological time series such as microarray & liquid/gas chromatography-based laboratory data (including mass spectrometry and ultraviolet diode arrays). Although alignment and difference detection have been extensively studied as separate problems in the signal processing and statistical pattern recognition communities, to our knowledge, no existing model performs both tasks in a unified way. Single class alignment algorithms attempt to align a set of time series all together, assuming that variability across different time series is attributable purely to noise. In many real-world situations, however, we have time series from multiple classes (categories) and our prior belief is that there is both substantial shared structure between the class distributions and, simultaneously, systematic (although often rare) differences between them. While in some circumstances (if differences are small and infrequent), single class alignment can be applied to multi-class data, it is much more desirable to have a model which performs true multi-class alignment in a principled way, allowing for more refined and accurate modeling of the data. In this paper, we introduce a novel hierarchical Bayesian model which simultaneously solves the multi-class alignment and difference detection tasks in a unified manner, as illustrated in Figure 1. The single-class alignment shown in this figure coerces the feature in region A for class 1 to be inappropriately collapsed in time, and the overall width of the main broad peak in class 2 to be inappropriately narrowed. In contrast, our multi-class model handles these features correctly. Furthermore, because our algorithm does inference for a fully probabilistic model, we are able to obtain quantitative measures of the posterior uncertainty in our results, which, unlike the point estimates produced by most current approaches, allow us to assess our relative confidence in differences learned by the model. Our basic setup for multi-class alignment assumes the class labels are known for each time series, as is the case for most difference detection problems. However, as we discuss at the end of the paper, our model can be extended to the completely unsupervised case. 3

We present a model for predicting HLA class I restricted CTL epitopes. In contrast to almost all other work in this area, we train a single model on epitopes from all HLA alleles and supertypes, yet retain the ability to make epitope predictions for specific HLA alleles. We are therefore able to leverage data across all HLA alleles and/or their supertypes, automatically learning what information should be shared and also how to combine allele-specific, supertype-specific, and global information in a principled way. We show that this leveraging can improve prediction of epitopes having HLA alleles with known supertypes, and dramatically increases our ability to predict epitopes having alleles which do not fall into any of the known supertypes. Our model, which is based on logistic regression, is simple to implement and understand, is solved by finding a single global maximum, and is more accurate (to our knowledge) than any other model. These keywords were added by machine and not by the authors. This process is experimental and the keywords may be updated as the learning algorithm improves.

Multiple realizations of continuous-valued time series from a stochastic process often contain systematic variations in rate and amplitude. To leverage the information contained in such noisy replicate sets, we need to align them in an appropriate way (for example, to allow the data to be properly combined by adaptive averaging). We present the Continuous Profile Model (CPM), a generative model in which each observed time series is a non-uniformly subsampled version of a single latent trace, to which local rescaling and additive noise are applied. After unsupervised training, the learned trace represents a canonical, high resolution fusion of all the replicates. As well, an alignment in time and scale of each observation to this trace can be found by inference in the model. We apply CPM to successfully align speech signals from multiple speakers and sets of Liquid Chromatography-Mass Spectrometry proteomic data. 1 A Profile Model for Continuous Data When observing multiple time series generated by a noisy, stochastic process, large sys- tematic sources of variability are often present. For example, within a set of nominally replicate time series, the time axes can be variously shifted, compressed and expanded, in complex, non-linear ways. Additionally, in some circumstances, the scale of the mea- sured data can vary systematically from one replicate to the next, and even within a given replicate. We propose a Continuous Profile Model (CPM) for simultaneously analyzing a set of such time series. In this model, each time series is generated as a noisy transformation of a single latent trace. The latent trace is an underlying, noiseless representation of the set of replicated, observable time series. Output time series are generated from this model by moving through a sequence of hidden states in a Markovian manner and emitting an observable value at each step, as in an HMM. Each hidden state corresponds to a particular location in the latent trace, and the emitted value from the state depends on the value of the latent trace at that position. To account for changes in the amplitude of the signals across and within replicates, the latent time states are augmented by a set of scale states, which control how the emission signal will be scaled relative to the value of the latent trace. During training, the latent trace is learned, as well as the transition probabilities controlling the Markovian evolution of the scale and time states and the overall noise level of the observed data. After training, the latent trace learned by the model represents a higher resolution 'fusion' of the experimental replicates. Figure 1 illustrate the model in action. Unaligned, Linear Warp Alignment and CPM Alignment 40 30 20 Amplitude 10 0 50 40 30 20 Amplitude 10 0 30 20 Amplitude 10 0 Time a) b) Figure 1: a) Top: ten replicated speech energy signals as described in Section 4), Middle: same signals, aligned using a linear warp with an offset, Bottom: aligned with CPM (the learned latent trace is also shown in cyan). b) Speech waveforms corresponding to energy signals in a), Top: unaligned originals, Bottom: aligned using CPM. 2 Defining the Continuous Profile Model (CPM) The CPM is generative model for a set of K time series, xk = (xk1, xk2, ..., xk ). The N k temporal sampling rate within each xk need not be uniform, nor must it be the same across the different xk. Constraints on the variability of the sampling rate are discussed at the end of this section. For notational convenience, we henceforth assume N k = N for all k, but this is not a requirement of the model. The CPM is set up as follows: We assume that there is a latent trace, z = (z1, z2, ..., zM ), a canonical representation of the set of noisy input replicate time series. Any given observed time series in the set is modeled as a non-uniformly subsampled version of the latent trace to which local scale transformations have been applied. Ideally, M would be infinite, or at least very large relative to N so that any experimental data could be mapped precisely to the correct underlying trace point. Aside from the computational impracticalities this would pose, great care to avoid overfitting would have to be taken. Thus in practice, we have used M = (2 + )N (double the resolution, plus some slack on each end) in our experiments and found this to be sufficient with < 0.2. Because the resolution of the latent trace is higher than that of the observed time series, experimental time can be made effectively to speed up or slow down by advancing along the latent trace in larger or smaller jumps. The subsampling and local scaling used during the generation of each observed time se- ries are determined by a sequence of hidden state variables. Let the state sequence for observation k be k. Each state in the state sequence maps to a time state/scale state pair: k { k, k}. Time states belong to the integer set (1..M ); scale states belong to an i i i ordered set (1..Q). (In our experiments we have used Q=7, evenly spaced scales in logarithmic space). States, k, and observation values, xk, are related by the emission i i probability distribution: A (xk|z) p(xk|k, z, , uk) N (xk; z kuk, ) k , where i i i i k i i i is the noise level of the observed data, N (a; b, c) denotes a Gaussian probability density for a with mean b and standard deviation c. The uk are real-valued scale parameters, one per observed time series, that correct for any overall scale difference between time series k and the latent trace. To fully specify our model we also need to define the state transition probabilities. We define the transitions between time states and between scale states separately, so that T k p( i|i-1) = p(i|i-1)pk (i|i-1). The constraint that time must move i-1,i forward, cannot stand still, and that it can jump ahead no more than J time states is en- forced. (In our experiments we used J = 3.) As well, we only allow scale state transitions between neighbouring scale states so that the local scale cannot jump arbitrarily. These constraints keep the number of legal transitions to a tractable computational size and work well in practice. Each observed time series has its own time transition probability dis- tribution to account for experiment-specific patterns. Both the time and scale transition probability distributions are given by multinomials: dk 1, ifa-b=1 dk2, ifa-b=2 . pk( i = a|i-1 = b) = ..dk, ifa-b=J J 0, otherwise s0, ifD(a,b)=0 s p( 1, if D(a, b) = 1 i = a|i-1 = b) = s1, ifD(a,b)=-1 0, otherwise where D(a, b) = 1 means that a is one scale state larger than b, and D(a, b) = -1 means that a is one scale state smaller than b, and D(a, b) = 0 means that a = b. The distributions are constrained by: J dk = 1 and 2s i=1 i 1 + s0 = 1. J determines the maximum allowable instantaneous speedup of one portion of a time series relative to another portion, within the same series or across different series. However, the length of time for which any series can move so rapidly is constrained by the length of the latent trace; thus the maximum overall ratio in speeds achievable by the model between any two entire time series is given by min(J , M ). N After training, one may examine either the latent trace or the alignment of each observable time series to the latent trace. Such alignments can be achieved by several methods, in- cluding use of the Viterbi algorithm to find the highest likelihood path through the hidden states [1], or sampling from the posterior over hidden state sequences. We found Viterbi alignments to work well in the experiments below; samples from the posterior looked quite similar. 3 Training with the Expectation-Maximization (EM) Algorithm As with HMMs, training with the EM algorithm (often referred to as Baum-Welch in the context of HMMs [1]), is a natural choice. In our model the E-Step is computed exactly using the Forward-Backward algorithm [1], which provides the posterior probability over states for each time point of every observed time series, k(i) p( s i = s|x) and also the pairwise state posteriors, s,t(i) p(i-1 = s, i = t|xk). The algorithm is modified only in that the emission probabilities depend on the latent trace as described in Section 2. The M-Step consists of a series of analytical updates to the various parameters as detailed below. Given the latent trace (and the emission and state transition probabilities), the complete log likelihood of K observed time series, xk, is given by Lp L + P. L is the likelihood term arising in a (conditional) HMM model, and can be obtained from the Forward-Backward algorithm. It is composed of the emission and state transition terms. P is the log prior (or penalty term), regularizing various aspects of the model parameters as explained below. These two terms are: K N N L log p(1) + log A (xk|z) + log T k (1) i i i-1,i k=1 i=1 i=2 -1 K P - (zj+1 - zj)2 + log D(dk|{k}) + log D(s }), v v v |{v (2) j=1 k=1 where p(1) are priors over the initial states. The first term in Equation 2 is a smoothing penalty on the latent trace, with controlling the amount of smoothing. k v and v are Dirichlet hyperprior parameters for the time and scale state transition probability distribu- tions respectively. These ensure that all non-zero transition probabilities remain non-zero. For the time state transitions, v {1, J } and kv corresponds to the pseudo-count data for the parameters d1, d2 . . . dJ . For the scale state transitions, v {0, 1} and k v corresponds to the pseudo-count data for the parameters s0 and s1. Letting S be the total number of possible states, that is, the number of elements in the cross-product of possible time states and possible scale states, the expected complete log likelihood is: K S K S N <Lp>=P + k(1) log T k + k(i) log A |z) + . . . s 0,s s s(xk i k=1 s=1 k=1 s=1 i=1 K S S N . . . + k (i) log T k s,s s,s k=1 s=1 s =1 i=2 using the notation T k 0 p( (i) (i) ,s 1 = s), and where k s and k are the posteriors over s,s states as defined above. Taking derivatives of this quantity with respect to each of the parameters and finding the critical points provides us with the M-Step update equations. In updating the latent trace z we obtain a system of M simultaneous equations, for j = 1..M : K N <Lp> (xk - z = 0 = k(i) i j uk s) - (4z z s suk j - 2zj-1 - 2zj+1) j 2 k=1 {s|s=j} i=1 For the cases j = 1, N , the terms 2zj-1 and 2zj+1, respectively, drop out. Considering all such equations we obtain a system of M equations in M unknowns. Each equation depends only linearly on three variables from the latent trace. Thus the solution is easily obtained numerically by solving a tridiagonal linear system. Analytic updates for 2 and uk are given by: S N k(i)(xk - z uk S z N k(i)xk 2 = s=1 i=1 s i s s)2 , uk = s=1 s s i=1 s i N S (z k(i) s=1 s s)2 N i=1 s Lastly, updates for the scale and state transition probability distributions are given by: k + S N k (i) v s=1 {s | - i=2 s,s dk = s s =v} v J k + J S N k (i) j=1 j j=1 s=1 {s | - i=2 s,s s s =j} + K S N k (i) j k=1 s=1 {s H(s,v)} i=2 s,s sv = 1 + K S N k (i) j=0 j k=1 s=1 {s H(s,1),H(s,0)} i=2 s,s where H(s, j) {s |s is exactly j scale states away from s}. Note that we do not normalize the Dirichlets, and omit the traditional minus one in the exponent: D(dk|{k}) = J (dk)kv }) = 1 (s v v and D(s v=1 v v |{v v=0 v )v . The M-Step updates uk, , and z are coupled. Thus we arbitrarily pick an order to update them and as one is updated, its new values are used in the updates for the coupled parameter updates that follow it. In our experiments we updated in the following order: , z, uk. The other two parameters, dkv and sv, are completely decoupled. 4 Experiments with Laboratory and Speech Data We have applied the CPM model to analyze several Liquid Chromatography - Mass Spec- trometry (LC-MS) data sets from an experimental biology laboratory. Mass spectrometry technology is currently being developed to advance the field of proteomics [2, 3]. A mass spectrometer takes a sample as input, for example, human blood serum, and produces a measure of the abundance of molecules that have particular mass/charge ratios. In pro- teomics the molecules in question are small protein fragments. From the pattern of abun- dance values one can hope to infer which proteins are present and in what quantity. For protein mixtures that are very complex, such as blood serum, a sample preparation step is used to physically separate parts of the sample on the basis of some property of the molecules, for example, hydrophobicity. This separation spreads out the parts over time so that at each unique time point a less complex mixture is fed into the mass spectrometer. The result is a two-dimensional time series spectrum with mass/charge on one axis and time of input to the mass spectrometer on the other. In our experiments we collapsed the data at each time point to one dimension by summing together abundance values over all mass/charge values. This one-dimensional data is referred to as the Total Ion Count (TIC). We discuss alternatives to this in the last section. After alignment of the TICs, we assessed the alignment of the LC-MS data by looking at both the TIC alignments, and also the cor- responding two-dimensional alignments of the non-collapsed data, which is where the true information lies. The first data set was a set of 13 replicates, each using protein extracted from lysed E. coli cells. Proteins were digested and subjected to capillary-scale LC-MS coupled on-line to an ion trap mass spectrometer. First we trained the model with no smoothing (i.e., = 0) on the 13 replicates. This provided nice alignments when viewed in both the TIC space and the full two-dimensional space. Next we used leave-one-out cross-validation on six of the replicates in order to choose a suitable value for . Because the uk and dkv are time series specific, we ran a restricted EM on the hold-out case to learn these parameters, holding the other parameters fixed at the values found from learning on the training set. Sixteen values of over five orders of magnitude, and also zero, were used. Note that we did not include the regularization likelihood term in the calculations of hold-out likelihood. One of the non-zero values was found to be optimal (statistically significant at a p=0.05 level using a paired sample t-test to compare it to no smoothing). Visually, there did not appear to be a difference between no regularization and the optimal value of , in either the TIC space or the full two-dimensional space. Figure 2 shows the alignments applied to the TICs and also the two-dimensional data, using the optimal value of . 8 x 10 Unaligned and Aligned Time Series 8 x 10 Replicate 5 10 9 Latent Trace Original Time Series Aligned Experimental Time Series 8 8 7 6 6 Amplitude 4 5 2 4 Amplitude 8 3 0 x 10 2 6 1 0 100 200 300 400 500 600 700 800 Residual 4 3 Time Jump From Previous State 2 2 1 Latent Space Amplitude Scale States 0 100 200 300 400 500 600 700 800 200 400 600 800 Time Latent Time a) b) c) d) Figure 2: Figure 2: a) Top: 13 Replicate pre-processed TICs as described in Section 4), Bottom: same as top, but aligned with CPM (the learned latent trace is also shown). b) The fifth TIC replicate aligned to the learned latent trace (inset shows the original, unaligned). Below are three strips showing, from top-to-bottom, i) the error residual, ii) the number of time states moved between every two states in the Viterbi alignment, and iii) the local scaling applied at each point in the alignment. c) A portion of the two-dimensional LC-MS data from replicates two (in red) and four (in green). d) Same as c), but after alignment (the same one dimensional alignment was applied to every Mass/Charge value). Marker lines labeled A to F show how time in c) was mapped to latent time using the Viterbi alignment. We also trained our model on five different sets of LC-MS data, each consisting of human blood serum. We used no smoothing and found the results visually similar in quality to the first data set. To ensure convergence to a good local optimum and to speed up training, we pre-processed the LC-MS data set by coarsely aligning and scaling each time series as follows: We 1) translated each time series so that the center of mass of each time series was aligned to the median center of mass over all time series, 2) scaled the abundance values such that the sum of abundance values in each time series was equal to the median sum of abundance values over all time series. We also used our model to align 10 speech signals, each an utterance of the same sentence spoken by a different speaker. The short-time energy (using a 30ms Hanning window) was computed every 8ms for each utterance and the resulting vectors were used as the input to CPM for alignment. The smoothing parameter was set to zero. For comparison, we also performed a linear warping of time with an offset. (i.e. each signal was translated so as to start at the same time, and the length of each signal was stretched or compressed so as to each occupy the same amount of time). Figure 1 shows the successful alignment of the speech signals by CPM and also the (unsuccessful) linear warp. Audio for this exam- ple can be heard at http://www.cs.toronto.edu/~jenn/alignmentStudy, which also contains some supplemental figures for the paper. Initialization for EM training was performed as follows: was set to 15% of the difference between the maximum and minimum values of the first time series. The latent trace was initialized to be the first observed time series, with Gaussian, zero-mean noise added, with standard deviation equal to . This was then upsampled by a factor of two by repeating every value twice in a row. The additional slack at either end of the latent trace was set to be the minimum value seen in the given time series. The uk were each set to one and the multinomial scale and state transition probabilities were set to be uniform. 5 Related Algorithms and Models Our proposed CPM has many similarities to Input/Output HMMs (IOHMMs), also called Conditional HMMs [4]. IOHMMs extend standard HMMs [1] by conditioning the emission and transition probabilities on an observed input sequence. Each component of the output sequence corresponds to a particular component of the input. Training of an IOHMM is supervised -- a mapping from an observed input sequence to output target sequence is learned. Our CPM also requires input and thus is also a type of conditional HMM. However, the input is unobserved (but crucially it is shared between all replicates) and hence learning is unsupervised in the CPM model. One could also take the alternative view that the CPM is simply an HMM with an extra set of parameters, the latent trace, that affect the emission probabilities and which are learned by the model. The CPM is similar in spirit to Profile HMMs which have been used with great success for discrete, multiple sequence alignment, modeling of protein families and their con- served structures, gene finding [5], among others. Profile HMM are HMMs augmented by constrained-transition 'Delete' and 'Insert' states, with the former emitting no observa- tions. Multiple sequences are provided to the Profile HMM during training and a summary of their shared statistical properties is contained in the resulting model. The development of Profile HMMs has provided a robust, statistical framework for reasoning about sets of related discrete sequence data. We put forth the CPM as a continuous data, conditional analogue. Many algorithms currently used for aligning continuous time series data are variations of Dynamic Time Warping (DTW) [6], a dynamic programming based approach which origi- nated in the speech recognition community as a robust distance measure between two time series. DTW works on pairs of time series, aligning one time series to a specified reference time series. DTW does not take in to account systematic variations in the amplitude of the signal. Our CPM can be viewed as a rich and robust extension of DTW that can be applied to many time series in parallel and which automatically uncovers the underlying template of the data. 6 Discussion and Conclusion We have introduced a generative model for sets of continuous, time series data. By training this model one can leverage information contained in noisy, replicated experimental data, and obtain a single, superior resolution 'fusion' of the data. We demonstrated successful use of this model on real data, but note that it could be applied to a wide range of problems involving time signals, for example, alignment of gene expression time profiles, alignment of temporal physiological signals, alignment of motion capture data, to name but a few. Certain assumptions of the model presented here may be violated under different ex- perimental conditions. For example, the Gaussian emission probabilities treat errors in large amplitudes in the same absolute terms as in smaller amplitudes, whereas in real- ity, it may be that the error scales with signal amplitude. Similarly, the penalty term - -1(z j=1 j+1 - zj )2 does not scale with the amplitude; this might result in the model arbitrarily preferring a lower amplitude latent trace. (However, in practice, we did not find this to be a problem.) One immediate and straight-forward extension to the model would be to allow the data at each time point to be a multi-dimensional feature vector rather than a scalar value. This could easily be realized by allowing the emission probabilities to be multi-dimensional. In this way a richer set of information could be used: either the raw, multi-dimensional feature vector, or some transformation of the feature vectors, for example, Principal Components Analysis. The rest of the model would be unchanged and each feature vector would move as a coherent piece. However, it might also be useful to allow different dimensions of the feature vector to be aligned differently. For example, with the LC-MS data, this might mean allowing different mass/charge peptides to be aligned differently at each time point. However, in its full generality, such a task would be extremely computational intense. A perhaps more interesting extension is to allow the model to work with non-replicate data. For example, suppose one had a set of LC-MS experiments from a set of cancer patients, and also a set from normal persons. It would be desirable to align the whole set of time series and also to have the model tease out the differences between them. One approach is to consider the model to be semi-supervised - the model is told the class membership of each training example. Then each class is assigned its own latent trace, and a penalty is introduced for any disagreements between the latent traces. Care needs to be taken to ensure that the penalty plateaus after a certain amount of disagreement between latent trace points, so that parts of the latent trace which are truly different are able to whole-heartedly disagree. Assuming that the time resolution in the observed time series is sufficiently high, one might also want to encourage the amount of disagreement over time to be Markovian. That is, if the previous time point disagreed with the other latent traces, then the current point should be more likely to disagree.