Robust dose-response curve estimation applied to high content screening data analysis
© Nguyen et al.; licensee BioMed Central. 2014
Received: 14 July 2014
Accepted: 14 November 2014
Published: 10 December 2014
Background and method
Successfully automated sigmoidal curve fitting is highly challenging when applied to large data sets. In this paper, we describe a robust algorithm for fitting sigmoid dose-response curves by estimating four parameters (floor, window, shift, and slope), together with the detection of outliers. We propose two improvements over current methods for curve fitting. The first one is the detection of outliers which is performed during the initialization step with correspondent adjustments of the derivative and error estimation functions. The second aspect is the enhancement of the weighting quality of data points using mean calculation in Tukey’s biweight function.
Results and conclusion
Automatic curve fitting of 19,236 dose-response experiments shows that our proposed method outperforms the current fitting methods provided by MATLAB®;’s nlinfit nlinfit function and GraphPad’s Prism software.
In recent years, the need for automatic data analysis has been amplified by the use of high content screening (HCS)  techniques. HCS (also known as phenotypic screening) is a great tool to identify small molecules that alter the disease state of a cell based on measuring cellular phenotype. However, HCS always comes with an important caveat: there is little or no a priori information on the compound’s target. At the Institut Pasteur Korea (IP-K), we have applied high content screening using small interfering RNA (siRNA) , at the genome scale. More recently, we applied the siRNA knockdown strategy for each gene from the genome and studied their effect on a given molecule producing a clear response on a cellular assay, to look for synergestic effects of a drug and target. Performing large-scale curve fitting for the biological activity of a large compound library can lead to a great number of dose-response curves (DRCs) that require proper fitting. The overall process of fitting, rejecting outliers, and data point weighting of thousands to millions of curves can be a significant challenge. The quality of a dose-response curve fitting algorithm is a key element that should help in reducing false positive hits and increasing real compound hits.
Several solutions - have been proposed to perform curve fitting. One of the widely used nonlinear curve fitting algorithms was introduced by Levenberg-Marquardt (LM) ,. This method belongs to the gradient-descent family. However, due to the sensitivity of the method, the data quality and the initial guess, the curve fitting algorithm is unfortunately susceptible to becoming trapped in local minimums. To obtain proper curve fitting results, there is a need for automatic outlier handling, usage of predefined initial curves, and adaptive weighting for data points ,. All of these demands can, without doubt, be applicable to large sets of data. Nonlinear and non-iterative least square regression analysis was presented in  for robust logistic curve fitting with detection of possible outliers. This non-iterative algorithm was implemented in a microcomputer and assessed using different biological and medical data. A review of popular fitting models using linear and nonlinear regression is given in . The study serves as a practical guide to help researchers who are not statisticians understand statistical tools more clearly, especially dose-response curve fitting using nonlinear regression. This book also describes the robust fitting algorithm and the outlier detection mechanism used in the commercial software GraphPad Prism®;. In , an automatic best-of-fit estimation procedure is introduced based on the Akaike information criterion. This best-of-fit model guarantees not only the estimation of the parameters with smallest sum-of-squares errors but also good prediction of the model.
Simulated data for DRC estimation was used in  to examine the performance of the proposed Grid algorithm. The peculiarity of the Grid algorithm is that it visits all the points in a grid of four curve parameters and searches for the point with the optimum sum-of-squares error. A coarse-to-fine grid model and a threshold-based outlier detection mechanism were used to make the algorithm more efficient. This paper also provided Java-based software together with a sample dataset for academic use. However, the software is applicable only when the measurements are taken without replicates (one data point at each concentration). Furthermore, various popular computer software and code packages have been presented for DRC fitting such as the DRC package in R , the nlinfit nlinfit function in MATLAB®; , and the XLfit add-in for Microsoft®; Excel®; .
In this paper, robust fitting and automatic outlier detection based on Tukey’s biweight function are introduced. This method was developed to automate nonlinear fitting of thousands of DRCs performed in replicates, detect the outliers automatically, and initialize the fitting curves robustly.
Background and method
Basic computation of curve fitting
Evaluate χ(β ) and define a modest value for λ, i.e., λ=0.001;
If , increase λ by a factor of 10; else, decrease λ by a factor of 10 and update β ←β +δ β ;
Repeat steps 2 and 3 until χ(β ) converges, and the return β .
instead of (8) which is used without outlier detection. Herein, we proposed (10) to reduce the impact of outliers on the fitting results by considering absolute errors and sign-only derivatives. The key difference between (10) and (8) is that the derivative function: ψ(z)=s i g n(z) results in −1 (negative) or 1 (positive), which controls the gradient in the direction of having a higher number of negatives/positives, whereas ψ(z)=z judges the gradient based on the distance between the estimated and actual values. Therefore, (10) is able to disregard the points having a lower number of negatives/positives (see Figure 2, on the left), and (8) considers the points at far distances although these points are given low impact (see Figure 2, on the right). Levenberg-Marquardt and other conventional nonlinear curve fitting algorithms are based on derivative calculation, and the quality of their solutions notably depends on data quality (i.e., outliers) and the initial guess. To have good fitting, outliers and the initial guess have to be manually detected and defined. Accordingly, these conventional algorithms are very difficult to automate and be made to yield good solutions in thousands of DRCs. Based on (10), outliers can be effectively weighted and a robust initial guess is automatically determined at the beginning of the fitting process. The method of weighting data points is described in the next section.
Curve fitting with weighting function
Because all fitting parameters of the DRC are very important in understanding and assessing the effect of a chemical compound, it is essential to have a method that can estimate the curve in a robust way, i.e., coping with outliers and assigning weights to data points. Initially all data points are supposed to have equal weights. However, this idea does not hold in many practical occasions. Therefore, a least squares method tends to give unequal weighting to data points, e.g., points that are closer to the fitted curve would have higher weighting values. In standard weighting, minimizing the fitting error (sum-of-squares) of the absolute vertical distances is not appropriate: points having high response values tend to have large deviations from the curve and so they contribute more to the sum-of-squares value. This weighting makes sense when the scattering of data is Gaussian and the standard deviation among replicates is approximately the same at each concentration.
To overcome the situation in which data spreads differently at concentrations, various weighting techniques are considered, including relative weighting, Poisson weighting, and observed variability-based weighting . Relative weighting extends the idea of standard weighting by dividing the squared distance by the square of the corresponding response value Y; hence, the relative variability is consistent. Similarly, Poisson weighting and weighting by observed variability use different forms of dividing the response value Y. Indeed, minimizing the sum-of-squares might yield the best possible curve when all variations obey a Gaussian distribution (without considering how different the standard deviations at concentrations are). However, it is common for one data point to be far from the rest (caused by experimental mistakes); then, this point does not belong to the same Gaussian distribution as the remaining points and it contributes erroneous impact to the fitting. The Tukey biweight function  was introduced to reduce the effect of outliers. This weighting function considers large residuals and treats them with low weights, or even zero weights, so that they do not sway the fitting much. In this section, we present a modification of the Tukey function and apply it to our fitting.
Determine the distance from each data point to the curve, called the residual, r;
Calculate the weight of each point using (11) with applied;
Assign new values to data points based on their weights.
In addition, other robust weighting functions can be considered such as Andrews, bisquare, Cauchy, Fair, Huber, logistic, Talwar, and Welsch ; however, the Tukey biweight function is known for its reliability, and it is generally recommended for robust optimization . As we show further, the use of the Tukey biweight allows for improved DRC fitting results (see the Results section).
Robust univariate DRC estimation algorithm
Based on the curve obtained, calculate the Tukey biweights of data points using (11);
Based on the obtained weights, execute the LM algorithm in the basic computation section and consider the weights of the data points.
Repeat steps 2 and 3 for a predefined number of iterations or until convergence.
We have previously shown the ability to use the high content genome-wide silencing RNA (siRNA) screening approach on various cellular models ,. The recent combination of drugs (at different concentrations) and siRNA (approximately 20,000) lead us to look for an automatic method to characterize a large number of dose-response curves obtained in those experimental conditions. Using MATLAB®;, a curve fitting algorithm was implemented. We compared the fitting performance of our method to that of the nlinfit nlinfit function in MATLAB®; 2013a and the robust DRC fitting in GraphPad Prism®; 6.0.
In our experiments, we used eight different concentrations 0, 0.005, 0.01, 0.1, 0.4, 1, 5, and 20 μ M with five replicates at concentration 0 and three replicates at the other concentrations (26 data points in total). In this manuscript, the concentrations were plotted over the X-axis in log unit. The Y-axis shows the drug response, which was normalized in the range of 0 to 100. The initial values of the parameters floor, window, shift, and slope used in the MATLAB®; nlinfit nlinfit function were defined based on values of data points as m i n(Y), m a x(Y)−m i n(Y), a v e r a g e(X), and –0.6, respectively. By default, the algorithm uses bisquare bisquare (also known as Tukey biweight) as the robust weighting function. Indeed, MATLAB®; applied the Levenberg-Marquardt (LM) algorithm and iterative reweighted least squares  for robust estimation. Accordingly, the nlinfit nlinfit represents the case of using the traditional LM algorithm and Tukey biweight function where (8) and the median function are utilized. For our algorithm, we defined the initial DRC parameters in the first step (outlier detection) as in nlinfitnlinfit, and then those parameters were corrected using the outlier detection step and used for the consequent fitting steps. According to step 4 of our DRC estimation algorithm, the number of iterations was predefined as 50, and the error tolerance for convergence was 0.0001.
Averages and standard deviations of the normalized sum-of-squares errors calculated based on the fitting results of 19,236 curves
Number of curves
In summary, experimental comparisons show that our method (namely ‘Ours (mean)’ in the figure), which proposes automatic initialization of DRC parameters and modification of the Tukey biweight function (mean calculation), yields a satisfactory fitting of curves. It provides more accurate fitting than MATLAB®; nlinfitnlinfit, where automatic initialization is not available and the default Tukey biweight function (median calculation) is used, by more than 14% in processing 19,236 curves. We also demonstrated that the method applying the automatic initialization and the default Tukey function (namely ‘Ours (median)’) did not yield results as good as those of ‘Ours (mean)’. Moreover, the better performance of ‘Ours (median)’ than that of MATLAB®; nlinfit nlinfit implies that the automatic initialization of DRC parameters meaningfully improves the fitting process. Our method is superior to GraphPad Prism®; 6.0 by more than 1%. Although the error is slightly improved; we believe that, with the MATLAB®; implementation provided, our approach is easily automated and scalable to thousands of curves. It is able to process the entire genome data in less than two hours – approximately 360 milliseconds for a curve (with a computer configuration using an Intel Core i7 3.47 GHz CPU).
We provide two improvements for the problem of DRC fitting: 1) increasing the accuracy of the initialization of DRC parameters with the use of outlier detection, and 2) improving the method of weighting for noisy data in the Tukey biweight function. Our method is adapted to the analysis of thousands of DRCs or more with the use of automatic outlier detection and initialization of curves. By experimentally comparing the results of our method to those calculated by the nlinfit nlinfit function in MATLAB®; 2013a and the robust DRC fitting in GraphPad Prism®; 6.0, we found that the proposed approach yielded a superior estimation of curves to that of MATLAB®; and Prism®;.
This work was supported by the National Research Foundation of Korea (NRF) grant funded by the Korea government (MEST) (No. 2012-00011), Gyeonggi-do and KISTI. Kyungmin Song and Myungjoo Kang were supported by the National Research Foundation of Korea (NRF) grant funded by the Korea government (MEST) (2013-025173). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.
- Joseph MZ:Applications of high content screening in life science research. Comb Chem High Throughput Screen. 2009, 12 (9): 870-876. 10.2174/138620709789383277.View ArticleGoogle Scholar
- Siqueira-Neto JL, Moon SH, Jang JY, Yang GS, Lee CB, Moon HK, Chatelain E, Genovesio A, Cechetto J, Freitas-Junior LH:An image-based high-content screening assay for compounds targeting intracellular Leishmania donovani amastigotes in human macrophages. PLOS Neglect Trop D. 2012, 6 (6): e1671-10.1371/journal.pntd.0001671.View ArticleGoogle Scholar
- Genovesio A, Kwon YJ, Windisch MP, Kim NY, Choi SY, Kim HC, Jung SY, Mammano F, Perrin V, Boese AS, Casartelli N, Schwartz O, Nehrbass U, Emans N:Automated genome-wide visual profiling of cellular proteins involved in HIV infection. J Biomol Screen. 2011, 16 (9): 945-958. 10.1177/1087057111415521.View ArticlePubMedGoogle Scholar
- Motulsky H, Christopoulos A: Fitting Models to Biological Data Using Linear and Nonlinear Regression: a Practical Guide to Curve Fitting . 2004, Oxford University Press Inc., New York, USAGoogle Scholar
- Levenberg K:A method for the solution of certain problems in least squares. Quart Applied Math. 1944, 2: 164-168.Google Scholar
- Marquardt D:An algorithm for least-squares estimation of nonlinear parameters. SIAM J Applied Math. 1963, 11 (2): 431-441. 10.1137/0111030.View ArticleGoogle Scholar
- Ayiomamitis A:Logistic curve fitting and parameter estimation using nonlinear noniterative least-squares regression analysis. Comput Biomed Res. 1986, 19 (2): 142-150. 10.1016/0010-4809(86)90012-1.View ArticlePubMedGoogle Scholar
- Rey D: Automatic best of fit estimation of dose response curve. Konferenz der SAS-Anwender in Forschung und Entwicklung 2007. [http://www.uni-ulm.de/ksfe2007]
- Wang Y, Jadhav A, Southal N, Huang R, Nguyen DT:A Grid algorithm for high throughput fitting of dose-response curve data. Curr Chem Genomics. 2010, 4: 57-66. 10.2174/1875397301004010057.PubMed CentralView ArticlePubMedGoogle Scholar
- Hoaglin D, Mosteller F, Tukey J: Understanding Robust and Exploratory Data Analysis . 1983, John Wiley and Sons Inc., New York, USAGoogle Scholar
- Boyd Y, Vandenberghe L: Convex Optimization . 2009, Cambridge University Press, CambridgeGoogle Scholar
- Ritz C, Streibig JC:Bioassay analysis using R. J Stat Soft. 2005, 12 (5): 1-22.View ArticleGoogle Scholar
- MathWorks Matlab 2013: nlinfit – Nonlinear regression[http://www.mathworks.com/help/stats/nlinfit.html]
- IDBS 2013: XLfit – for curve fitting and data analysis[http://www.excelcurvefitting.com]
- Press WH, Teukolsky SA, Vetterling WT, Flannery BP: Numerical Recipes in C: the Art of Scientific Computing (3rd edition) . 2007, Cambridge University Press, New York, USAGoogle Scholar
- Holland PW, Welsch RE:Robust regression using iteratively reweighted least-squares. Comm Stat Theor Meth. 1977, A6: 813-827. 10.1080/03610927708827533.View ArticleGoogle Scholar
- Maronna R, Martin RD, Yohai V: Robust Statistics – Theory and Methods . 2006, Wiley, Chichester, EnglandView ArticleGoogle Scholar
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly credited. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.