Mixture models for analysis of melting temperature data

Background In addition to their use in detecting undesired real-time PCR products, melting temperatures are useful for detecting variations in the desired target sequences. Methodological improvements in recent years allow the generation of high-resolution melting-temperature (Tm) data. However, there is currently no convention on how to statistically analyze such high-resolution Tm data. Results Mixture model analysis was applied to Tm data. Models were selected based on Akaike's information criterion. Mixture model analysis correctly identified categories in Tm data obtained for known plasmid targets. Using simulated data, we investigated the number of observations required for model construction. The precision of the reported mixing proportions from data fitted to a preconstructed model was also evaluated. Conclusion Mixture model analysis of Tm data allows the minimum number of different sequences in a set of amplicons and their relative frequencies to be determined. This approach allows Tm data to be analyzed, classified, and compared in an unbiased manner.


Background
Real-time PCR or semiquantitative PCR is widely used to detect and quantify specific target sequences. The exponential amplification of a sequence is monitored in real time by fluorescence. Commonly, a nonspecific fluorescent dye is used, such as SYBR Green I or LCGreen, which only reports the presence of double-stranded DNA. These dyes do not distinguish sequences and can thus report the amplification of undesired targets. Undesired sequences are normally detected during a dissociation step after thermocycling is complete. During dissociation, the doublestranded PCR products melt into single strands, so fluorescence is diminished. A curve can be produced by plot-ting the loss of fluorescence against a gradual increase in temperature. The temperature at which the rate of signal loss is the greatest can be defined as the melting temperature (T m ) of the PCR product. Although the T m is sequence dependent, different sequences do not necessarily have different T m . However, the converse is true. The detection of different T m does imply the presence of different sequences. Therefore, by monitoring T m , we can distinguish different targets for one set of primers. This technique has been used for the detection of single-nucleotide polymorphisms [1], allelic discrimination [2], and strain typing of microorganisms [3][4][5]. We previously reported the use of T m analysis to detect the expression patterns of transcripts containing different members of the W family of human endogenous retrovirus (HERV) elements in vitro and in vivo [6,7].
The precision of the T m measurements determines the sensitivity with which different sequences can be distinguished. The instrument used to obtain the T m recordings is the principal factor limiting the amount of information that can be extracted from the data. We recently reported a method that allows improved resolution, reduced spatial bias, and automated data collection for T m detection in an ABI Prism 7000 Sequence Detection System (Applied Biosystems, Palo Alto, CA) [8]. Using a temperature indicator probe (T m probe) and an algorithm (GcTm) to interpolate more-precise T m measurements from multiple data points, the standard deviation of the measurement error (σ) of the T m recordings was improved from 0.19°C to 0.06°C [8].
However, there is no convention on how to analyze T m data to objectively distinguish sequences by T m . The need for such a tool becomes apparent when the T m data are: i) not easily stratified because of overlapping clusters of T m observations, and/or ii) if the number of different sequences and possible T m categories are unknown. In this report, we use mixture model analysis to construct a model for a particular set of primer targets, to classify T m data, and to calculate the mixing proportions of the amplicons within these categories. The mixture model technique allows T m analysis to be applied to any set of primers to determine the minimum number of T m categories (i.e., the number of different sequences detected) and the mixing proportions (frequency distributions) of the detected categories. Thus, mixture model analysis of T m data is an objective method with which more refined T m assays can be established.

Results
In a T m analysis using the T m probe and GcTm program, described previously [8], we demonstrated, using plasmids containing known sequences, that it was possible to distinguish some but not all sequences based on their T m . In the present report, we applied the mixture models and the ρ established in the previous publication [8] to determine the T m categories and mixing proportions of these data ( Figure 1). Akaike's information criterion (AIC), a measure of how well a model explains the data, with a penalty for the number of parameters estimated, determined that the T m of the four sequences were best represented by a three-category mixture model. This model precisely estimated the mixing proportions of the T m into the categories, attributing the correct number of T m recordings to each of the four sequences (where two of them shared a category). For an overview of the procedure for using mixture models to analyze T m data, see the Methods section.  We next assessed the performance of the mixture model analysis in constructing models for categories of T m with varying separations. Therefore, we generated simulated data points mimicking the T m of four sequences separated by multiples of σ. These data were used to identify the model that best explained the data according to AIC (see an example of the AIC plot in Figure 2) for a range of T m separations and numbers of data points ( Figure 3). Next, we evaluated the fit of the data points to preestablished models. For this purpose, we generated data points corresponding to a sample containing three of four possible T m represented in a model. We compared the mixing proportions reported by the mixture model analysis with the mixing proportions in which all four T m were present at equal frequencies. In Figure 4, the P values obtained from χ 2 analyses for various separations of the T m are plotted against the numbers of data points used. The P values for the χ 2 test drop rapidly with increasing sample numbers for any T m separation of more than 1 × σ. With smaller separations of the T m categories, the mixture model analysis is unable to reliably establish the differences in the mixing proportions.

Discussion
We report the application of mixture models to the analysis of high-resolution T m data. Whereas the plasmid T m data reported are sufficiently separated to be stratified manually, we use these data to demonstrate the principle that can be applied to analyze more complex T m data.
Mixture model analysis of T m data entails the construction of a model based on the T m data for a set of primers. With such a model established, it is possible to fit smaller subsets of data to calculate the mixing proportions of the T m categories of the model. This gives a proxy marker for the frequency distributions of different amplicon sequences in the analyzed data. This approach requires no prior knowledge of how many different amplicons are present and there is no limit to the number of different T m that can be distinguished. However, the T m analysis method with mixture models only reports the minimum number of different sequences required to explain the T m data because different sequences can have the same T m .
Mixture model analysis is a modern type of cluster analysis. The purpose of cluster analysis is to group data that have properties in common. When constructing the mixture model for a set of primers, the number of categories in the model that most appropriately explains the T m data is determined by AIC. Other information criteria exist, such as the Bayesian information criterion, but this penalizes free parameters more harshly than does the AIC.
By empirical testing with simulated data, we found that smaller separations of T m require exponentially larger numbers of data points to distinguish the correct number of categories in a mixture model. Insufficient numbers of observations yield an underestimation of the numbers of unique T m represented by the data, erring on the side of safety. In other words, with insufficient data, the number of unique sequences in the data is underestimated by the optimal model.  AIC score Categories in model by AIC. In other words, models constructed with mixture model analysis will consist of T m categories separated by more than 1 × σ.
Not all dissociation curves are easily defined by a single T m , as in the case of multiple domain transitions in longer sequences [9] (generally longer than those generated in real-time PCR assays) and for heterodimers. Using the GcTm approach to curve fitting and SYBR Green I chemistry, such melting profiles will be assigned a single T m value. Although some additional information is therefore lost, mixture model analysis still validly identifies clusters of T m and sequences. There is an established high-resolution amplicon melting analysis (usually denoted HRM) using LCGreen, primarily based on differences in the profiles of melting curves rather than on absolute T m [10]. Although this method is superior to mixture model analysis in identifying heterodimers, absolute T m values are required to identify homodimers. Recently, a method with sufficient resolution to distinguish base-pair neutral homozygotes was reported [11]. Mixture model analysis of T m can be used in all cases where the T m can be denoted as a single value, but primarily for homodimer discrimination.

Conclusion
In conclusion, the mixture model analysis of T m presented here allows the unbiased analysis of high-resolution T m data. This analysis is applicable to the identification of sequences in T m data regardless of the method by which the T m are acquired, provided the measurement error is known. Mixture models allow T m analyses to be performed on more complex and varied sequence targets than hitherto possible. Possible applications include typing microbial strains and their relative abundances in a population and the analysis of transcripts containing repetitive elements [3,4,6,12].

Finite mixture models
Mixture models are useful for describing complex populations with observed or unobserved heterogeneity. The Plot of the optimal model, as defined by AIC, for simulated data consisting of four T m with varying separations against the number of T m used to construct the model Let X be a random variable or vector taking values in sample space χ with the probability density function where 0 ≤ π i ≤ 1, i = 1, ..., k, π 1 + ... + π k = 1.
Such a model can arise if one is sampling from a heterogeneous population that can be decomposed into k distinct homogeneous subpopulations, called component populations. If these components have been "mixed" together, and we measure only the variable X without determining the particular components, then this model holds. We say that X has a finite mixture distribution and that g(·) is a finite mixture density function. The parameters π 1 ,..., π k are called mixing weights or mixing proportions, and each π i represents the proportion of the total population in the i-th component.
There is no requirement that the component densities should all belong to the same parametric family, but in this paper, we keep to the simplest case where f 1 (x),..., f k (x) have a common functional form but different parameters.
We apply the theory of finite mixture models to T m data consisting of normally distributed components in a mixture model, where each component has a standard deviation of σ°C. The finite mixture density function is then as follows: where ψ = (π 1 ,..., π k , μ 1 ,..., μ k , σ) T .
The likelihood function corresponding to the data (x 1 ,..., x n ) is as follows: The logarithm of the likelihood function is We attempt to find the particular ψ that maximizes the likelihood function. This maximization can be undertaken in the traditional way by differentiating L(ψ; ×) with respect to the components of ψ and equating the derivatives to zero to give the likelihood equation: Quite often, the log likelihood function cannot be maximized analytically, i.e., the likelihood equation has no explicit solutions. In such cases, it is possible to compute the maximum likelihood of ψ iteratively. To calculate maximum likelihood estimates, we use the expectation maximization (EM) method in combination with the Newton-Raphson algorithm. Iterations of the EM algorithm consist of two steps: the expectation step or the E-step and the maximization step or the M-step [13,14]. The Newton-Raphson algorithm for solving the likelihood equation approximates the gradient vector of the log likelihood function by a linear Taylor series expansion [15]. We use the Newton-Raphson algorithm in the Mstep of the EM method.
We developed an algorithm that allows the automated estimation, in parallel, of a finite number of normally distributed components. The number of components can be assessed by several different methods, although none of them is optimal. We chose the AIC [16,17]. AIC is a relative score between different models where the selection of Separation (x σ) p -value Number of T data points m the optimal model is made by considering the number of data points and categories and the separation of the T m categories. AIC is defined as -2L m + 2m, where L m is the maximized log likelihood and m is the number of parameters.

Acquisition of HERV-W gag T m
T m data were generated with GcTm, as previously described [8], on dissociation data obtained from the amplification of plasmids containing known HERV-W gag sequences.
Simulated T m data recordings and GcTm analysis were performed in MATLAB™ (The MathWorks) version 7.0.1.24704 with the Optimization Toolbox. Mixture model analysis was performed in R 2.6.0 [18] with the MIX software [19,20].

Overview of mixture model analysis of T m
A mixture model is constructed for a set of primers. The model should be constructed on a large enough sample of T m data to expect all possible sequences to be represented. The T m data are then stratified into small-interval groups and the frequency distributions of these arbitrary categories are used to construct and compare the mixture models. AIC is used to evaluate which model best explains the data, while a minimum number of different categories is used. Lower values of AIC indicate the preferred model, i.e., the one with the fewest parameters. Once a model is selected, T m data from different samples can be fitted to the model and the mixing proportions compared between samples. Differences between samples can be evaluated with χ 2 tests if a conservative stance is taken, depending on the separation between the T m categories and the numbers of data points.

Authors' contributions
CN conceived the study, tested and prepared the manuscript; FU developed the method and critically revised the manuscript; JT developed the method and prepared the manuscript; HK conceived the study and prepared the manuscript.