Dorothee Childs1,2, Karsten Bach3,4, Holger Franken2, Simon Anders5, Nils Kurzawa1, Marcus Bantscheff2, Mikhail M Savitski1, Wolfgang Huber6. 1. European Molecular Biology Laboratory (EMBL) Heidelberg, Meyerhofstraβe 1, 69117 Heidelberg, Germany. 2. Cellzome GmbH, GlaxoSmithKline, Meyerhofstraβe 1, 69117 Heidelberg, Germany. 3. Department of Pharmacology, University of Cambridge, CB2 1PD, Cambridge, UK. 4. Cancer Research UK Cambridge Cancer Centre, CB2 0RE, Cambridge, UK. 5. Center for Molecular Biology of Heidelberg University (ZMBH), Im Neuenheimer Feld 282, 69120 Heidelberg, Germany. 6. European Molecular Biology Laboratory (EMBL) Heidelberg, Meyerhofstraβe 1, 69117 Heidelberg, Germany wolfgang.huber@embl.org.
Abstract
Detecting the targets of drugs and other molecules in intact cellular contexts is a major objective in drug discovery and in biology more broadly. Thermal proteome profiling (TPP) pursues this aim at proteome-wide scale by inferring target engagement from its effects on temperature-dependent protein denaturation. However, a key challenge of TPP is the statistical analysis of the measured melting curves with controlled false discovery rates at high proteome coverage and detection power. We present nonparametric analysis of response curves (NPARC), a statistical method for TPP based on functional data analysis and nonlinear regression. We evaluate NPARC on five independent TPP data sets and observe that it is able to detect subtle changes in any region of the melting curves, reliably detects the known targets, and outperforms a melting point-centric, single-parameter fitting approach in terms of specificity and sensitivity. NPARC can be combined with established analysis of variance (ANOVA) statistics and enables flexible, factorial experimental designs and replication levels. An open source software implementation of NPARC is provided.
Detecting the targets of drugs and other molecules in intact cellular contexts is a major objective in drug discovery and in biology more broadly. Thermal proteome profiling (TPP) pursues this aim at proteome-wide scale by inferring target engagement from its effects on temperature-dependent protein denaturation. However, a key challenge of TPP is the statistical analysis of the measured melting curves with controlled false discovery rates at high proteome coverage and detection power. We present nonparametric analysis of response curves (NPARC), a statistical method for TPP based on functional data analysis and nonlinear regression. We evaluate NPARC on five independent TPP data sets and observe that it is able to detect subtle changes in any region of the melting curves, reliably detects the known targets, and outperforms a melting point-centric, single-parameter fitting approach in terms of specificity and sensitivity. NPARC can be combined with established analysis of variance (ANOVA) statistics and enables flexible, factorial experimental designs and replication levels. An open source software implementation of NPARC is provided.
Keywords:
Drug targets; algorithms; biostatistics; functional data analysis; mathematical modeling; tandem mass spectrometry; thermal proteome profiling
Determining the cellular interaction partners of drugs and other small molecules remains a key challenge (1, 2, 3, 4). In drug research, better assays to detect targets (and off-targets) would provide valuable information on drugs' mechanisms of action, reveal potential reasons for side effects, and elucidate avenues for drug repurposing. More broadly, in cell biology basic research, the dynamical landscape of binding partners of metabolites, messengers or chemical probes contains much uncharted territory. Thermal proteome profiling (TPP) addresses these needs by screening for protein targets of drugs or small molecules in living cells on a proteome-wide scale (5, 6). TPP combines multiplexed quantitative mass spectrometry with the cellular thermal shift assay (CETSA) (7), which identifies binding events from shifts in protein thermostability (see supplemental Fig. S1 for a detailed explanation). A typical TPP experiment generates temperature dependent abundance measurements for a large part of the cellular proteome. Drug binding proteins can then be inferred by comparing the melting curves of proteins between samples treated with drug and vehicle (negative control without drug).Applications of TPP successfully identified previously unknown protein-ligand interactions (5), protein complexes (8) and downstream effects of drugs in signaling networks (6, 9, 10, 11) in human cells. Recently it has also been extended to study drug resistance in bacteria (12) and targets of antimalarial drugs in plasmodium (13). There is urgent interest in further advancing its component technologies, including experimental and computational aspects, in order to maximize its biological discovery potential (14, 15, 16, 17, 18, 19).The central computational task in TPP data analysis is the comparison of the temperature dependent abundance measurements—which can be visualized as melting curves— for each protein with and without (or with various concentrations of) drug. The aim is to detect changes in thermostability from statistically significant changes in the melting curves.A naïve approach is to summarize each curve into a single parameter, such as the melting point (Tm), which is defined as the temperature of half-maximum relative abundance (horizontal line in Fig. 1A). Its value is estimated by fitting a parametric model separately for the control and treatment conditions and comparing the estimates. Statistical significance is assessed using replicates and hypothesis testing, such as a t- or z-test. Although the approach has delivered valid and important results (5, 6, 20, 21), we will see in the following that it tends to lead to needlessly high rates of false negatives. There are three main reasons for that: first, drug-induced effects on thermostability do not always imply significant shifts in Tm (Fig. 1B–1C). Second, the true Tm of a protein can lie outside of the measured temperature range, which impairs its estimation (Fig. 1D). Both scenarios can result in important targets being missed in the analysis (Fig. 1E). The third reason is a more subtle statistical one: hypothesis tests using only the point estimates of Tm do not consider goodness-of-fit of the parametric model or the confidence range of the estimates. Thus, important information is ignored, which statistically leads to loss of power.
Fig. 1.
TPP data analysis challenges.
A–D, Examples for protein melting curves with and without drug (see color keys). In each case, ten temperatures were assayed, and two experimental replicates were made per condition, indicated by circle and triangle symbols. Fits of the sigmoid model (Eq. (3)) to both replicates jointly are shown by smooth lines. A, For serine/threonine protein kinase 4 (STK4), the binding of staurosporine is reflected by a marked shift between the curves. The fitted values for the melting points (Tm) are shown. B, For Bruton's tyrosine kinase (BTK), there is a small but reproducible shift between the curves. C, Protein kinase C beta (PRKCB) is destabilized by staurosporine; the effect occurs mainly at lower temperatures. D, NAD(P)H quinone dehydrogenase 2 (NQO2) is strongly stabilized by staurosporine. Although in each case, the effects of drug binding are clearly reproducible between replicates, the Tm-based approach of (6) only detects (A) and misses (B–D). In the case of (B) and (C), the fitted Tm are too similar, so that the statistical test does not assess the difference as significant. In the case of (D), no reasonable estimate for Tm in the staurosporine treated condition can be obtained, as it would lie outside the measured temperature range, and the protein is discarded from the analysis. In contrast, NPARC, the method proposed in this article, detects all four cases. E, The fraction of proteins in each of the data sets of Table I that is missed by the Tm-based approach because of failure to estimate Tm or to meet the goodness-of-fit criterion (Table III).
TPP data analysis challenges.
A–D, Examples for protein melting curves with and without drug (see color keys). In each case, ten temperatures were assayed, and two experimental replicates were made per condition, indicated by circle and triangle symbols. Fits of the sigmoid model (Eq. (3)) to both replicates jointly are shown by smooth lines. A, For serine/threonine protein kinase 4 (STK4), the binding of staurosporine is reflected by a marked shift between the curves. The fitted values for the melting points (Tm) are shown. B, For Bruton's tyrosine kinase (BTK), there is a small but reproducible shift between the curves. C, Protein kinase C beta (PRKCB) is destabilized by staurosporine; the effect occurs mainly at lower temperatures. D, NAD(P)H quinone dehydrogenase 2 (NQO2) is strongly stabilized by staurosporine. Although in each case, the effects of drug binding are clearly reproducible between replicates, the Tm-based approach of (6) only detects (A) and misses (B–D). In the case of (B) and (C), the fitted Tm are too similar, so that the statistical test does not assess the difference as significant. In the case of (D), no reasonable estimate for Tm in the staurosporine treated condition can be obtained, as it would lie outside the measured temperature range, and the protein is discarded from the analysis. In contrast, NPARC, the method proposed in this article, detects all four cases. E, The fraction of proteins in each of the data sets of Table I that is missed by the Tm-based approach because of failure to estimate Tm or to meet the goodness-of-fit criterion (Table III).
Table I
Datasets and sample sizes
Data set
Treatment
Concentration
Buffer
Cell line
Intact cells or lysate
Proteins
Reference
ATP data
MgATP
2 μm
PBS
K562
Lysate
4177
[9]
Dasatinib 0.5 μm data
Dasatinib
0.5 μm
PBS
K562
Intact Cells
4625
[5]
Dasatinib 5 μm data
Dasatinib
5 μm
PBS
K562
Intact Cells
4154
[5]
Panobinostat data
Panobinostat
1 μm
PBS
K562
Intact Cells
3649
[6]
Staurosporine data
Staurosporine
20 μm
PBS
K562
Lysate
4505
[5]
Table III
A priori filters applied in the original TPP analysis workflow [6] to select proteins for hypothesis testing
Rule number
Rule
1
Both fitted curves for the vehicle and compound treated condition have a coefficient of determination R2 > 0.8, where R2: = 1 − ∑t(xt−μc(t))2∑t(xt−x¯)2, with μc(t) being the model prediction for condition c at temperature t, x̄ being the mean of all measurements for the protein within a particular condition and replicate, and the summation extending over all temperatures.
2
The two curves fitted to the two replicates of the vehicle conditions have a plateau f∞ < 0.3.
3
In each biological replicate, the steepest slope of the melting curve in the vehicle and treatment condition needs to be < −0.06 °C−1.
Here, we propose an alternative approach that compares whole curves instead of summary parameters and does not rely on Tm estimation. The method, nonparametric analysis of response curves (NPARC), is based on a branch of statistical data analysis that works on continuous functions rather than individual numbers, termed functional data analysis (22). It considers the measured melting curves as samples from an underlying stochastic process with a smooth mean function—which can be modeled parametrically or nonparametrically (23)—and constructs its hypothesis tests directly on these samples. NPARC's F-statistic uses a more flexible model that makes fewer assumptions on the data than Tm-estimation, is computationally more stable, and it directly uses the information from replicates. Consequently, reliable estimates of the null distribution of this statistic can be obtained, it shows higher sensitivity for small but reproducible effects, and failures because of model misspecification or outliers are reduced. This increases proteome coverage, which can make the difference between missing or detecting an important drug target.We demonstrate NPARC on the five published data sets introduced in Table I. We also compare its results to those of the Tm-based method used by (6). Three of the experiments used the cancer drugs panobinostat or dasatinib in different concentrations, one investigated the effects of the high-affinity, ATP-competitive pan-kinase inhibitor staurosporine, one the cellular metabolite ATP. Although the cancer drugs interact with limited sets of proteins, the two other compounds are promiscuous binders and affect the thermostability of a large fraction of the cellular proteome.
EXPERIMENTAL PROCEDURES
Data Sets and Preprocessing
Five TPP data sets (Table I) were obtained from the supplements of the respective publications. Each data set contained relative abundance measurements per protein and temperature which had been scaled to the value measured at 37 °C (the lowest of the ten temperatures assayed) and subjected to the global normalization procedure described by Savitski et al. (5). Only proteins quantified with at least one unique peptide in each of two replicates of the vehicle and compound treated conditions were included in the analysis; the resulting proteome coverages are listed in Table I.
Curation of Lists of Expected Targets
Lists of expected protein targets for the pan-kinase inhibitor staurosporine and ATP were obtained from Gene Ontology Consortium annotations via the Bioconductor annotation packages AnnotationDbi (version 1.36.2), org.Hs.eg.db (version 3.4.0) and GO.db (version 3.4.0). Terms and numbers of annotated proteins are shown in Table II.
Table II
Expected targets per dataset
Data set
Gene Ontology term
Number of proteins in data set with this term
ATP
ATP-binding
558
Staurosporine
protein kinase activity
187
Mathematical Model
NPARC is based on fitting two competing models to the data, a null model and an alternative model. The null model states that the relative protein abundance at temperature t (given in °C) is explained by a single smooth function μ0(t) irrespective of the treatment condition (Fig. 2A). The alternative model posits two condition-specific functions: μ(t) for the treatment condition and μ(t) for the vehicle condition (Fig. 2B). Deviations between observed data and fitted model are referred to as residuals, and the sum of squared residuals (RSS) serves as an indicator of each model's goodness-of-fit. We then compute
where x is the measured value at temperature t for experimental replicate i and condition c ∈ {V,T}, and the summations extend over all temperatures, replicates, and conditions.
Fig. 2.
Principles of NPARC, illustrated for protein STK4 under staurosporine treatment.
A, Fit of the null model, i.e. no treatment effect (black line). The goodness-of-fit is quantified by RSS0, the sum of squared residuals (dashed lines). As in Fig. 1, the triangle and circle symbols indicate the experimental replicates. B, Fit of the alternative model, with separate curves for the treated (orange) and the vehicle condition (gray). Because of the higher flexibility of the model, the sum of squared residuals RSS1 is always less than or equal to RSS0. C, The question whether the improvement in the goodness-of-fit, i.e. the difference RSS0 − RSS1, is strong enough to reject the null hypothesis can be addressed with the variant of the F-test described in the main text. Each point in the plot corresponds to a different protein. The highlighted example STK4 has a large F-statistic and a small p value.
Principles of NPARC, illustrated for protein STK4 under staurosporine treatment.
A, Fit of the null model, i.e. no treatment effect (black line). The goodness-of-fit is quantified by RSS0, the sum of squared residuals (dashed lines). As in Fig. 1, the triangle and circle symbols indicate the experimental replicates. B, Fit of the alternative model, with separate curves for the treated (orange) and the vehicle condition (gray). Because of the higher flexibility of the model, the sum of squared residuals RSS1 is always less than or equal to RSS0. C, The question whether the improvement in the goodness-of-fit, i.e. the difference RSS0 − RSS1, is strong enough to reject the null hypothesis can be addressed with the variant of the F-test described in the main text. Each point in the plot corresponds to a different protein. The highlighted example STK4 has a large F-statistic and a small p value.
Choice of the Mean Function
The mean functions μ0(t), μ(t) and μ(t) are each chosen from the same space of smooth functions f : ℝ+ → [0,1] spanned by the three parameters a, b ∈ ℝ, f∞ ∈ [0,1] and the prescriptionThe shape of these functions is sigmoid, and the functional form (Eq. 3) can be motivated by simplifying protein thermodynamics considerations (5). The mean functions and the RSS values are computed separately from the data for each protein. In order not to overburden the notation, we omit the protein indices.
Hypothesis Test Statistic
To discriminate between null and alternative models, we compute the F-statistic
with d2/d1 > 0 defined as below. F quantifies the relative reduction in residuals from null to alternative model. Although F is by definition always positive, it will be small for proteins not affected by the treatment, whereas a high value of F indicates a reproducible change in thermostability.
Null Distribution
To compute a p value from a value of the F-statistic (Eq. 4), we need its null distribution, i.e. its statistical distribution if the data generating process is described by a common mean function μ). If the residuals were independent and identically normal distributed, this distribution would be given by an analytical formula, namely that of the F(d1, d2)-distribution with parameters d1, d2 > 0, and these parameters—sometimes called degrees of freedom—would be explicitly given from the number of measurements and number of model parameters that go into the computation of RSS0 and RSS1. In practice, this is not the case, because the residuals are heteroscedastic (i.e. have different variances at different temperatures) and correlated. However, the family of F(d1, d2)-distributions is quite flexible, and we can approximate the distribution of the F-statistic (Eq. 4) on data occurring in practice with an F(d1, d2)-distribution with different “effective degrees of freedom” d1, d2. To this end, we separately approximate the numerator and denominator of F as
and use the fact that the ratio of two χ-distributed random variables χ1), χ2) has an F(d1, d2)-distribution (24). The scale parameter σ02 and the effective degrees of freedom d1 and d2 are estimated from the empirical distributions—across proteins—of RSS0 and RSS1. Thus, we assume that σ02, d1, d2 are the same for all proteins. We estimate σ02 from the moments of RSS0 − RSS1 as
where mean and variance are computed across proteins on the observed values of RSS0 − RSS1 (see Supplementary Methods for details). Then, d1 and d2 are obtained by numerical optimization of the likelihoods for models (Eq. 5) and (Eq. 6) using the fitdistr function of the R package MASS (25).
p Values
For each protein, a p value is computed from its F-statistic and the cumulative F-distribution with parameters d1, d2 as described above. The multiset of p values across all proteins is corrected for multiple testing with the method of Benjamini and Hochberg (26). The outcome of such an analysis is exemplarily shown in Fig. 2C.
RESULTS
Application to Panobinostat
We assessed the ability of NPARC to detect drug targets on a data set on panobinostat (Table I). Panobinostat is a broad-spectrum histone deacetylase (HDAC) inhibitor known to interact with HDAC1, HDAC2, HDAC6, HDAC8, HDAC10, and tetratricopeptide repeat protein 38 (TTC38) (6).Out of 3649 proteins reproducibly quantified across both biological replicates in both treatment conditions, NPARC yielded 16 proteins with Benjamini-Hochberg adjusted p values ≤ 0.01. They contained the expected HDAC targets (Fig. 3A–3E) as well as TTC38, the histone proteins H2AFV or H2AFZ (the two variants could not be distinguished by mass spectrometry), and zinc finger FYVE domain-containing protein 28 (ZFYVE28) (Fig. 3F–3H). These proteins were previously identified as direct or indirect targets of panobinostat (6, 11). In addition, eight more proteins were detected for which no direct or indirect interactions with panobinostat have been described (supplemental Fig. S3). They reached statistical significance because they either showed effect sizes comparable to known panobinostat targets, or more subtle but highly reproducible changes in a similar strength to those already described for dasatinib target BTK (Fig. 1B). We reanalyzed the more recent 2D-TPP data set of short-term (15 min) panobinostat-treatment of HepG2 cells (11) for these proteins. All of them were identified and quantified at sufficient peptide coverage, but none of them showed stabilization. We thus conclude that the additionally found proteins are likely not direct binders of panobinostat, but rather indirect effects, like altered protein-protein interactions or post-translational modifications. The longer (5 h) incubation time of the assay used to generate the panobinostat data set in Table I makes it more sensitive to such effects.
Fig. 3.
Direct and indirect targets of the HDAC inhibitor panobinostat detected by NPARC (FDR ≤ 0.01 according to the method of Benjamini and Hochberg (
A–E, Data and curve fits for five HDACs that show significant shifts in their thermostability. HDAC1 and HDAC2 are not detected by the Tm-based approach of (6), because the higher variance between the replicates of the panobinostat-treated condition leads to them being eliminated by the filter heuristics of that method (Table III). In contrast, NPARC naturally takes the variance into account in the computation of the F-statistic and does not require such filtering steps. F–H, Data and curve fits for known non-HDAC targets.
Direct and indirect targets of the HDAC inhibitor panobinostat detected by NPARC (FDR ≤ 0.01 according to the method of Benjamini and Hochberg (
A–E, Data and curve fits for five HDACs that show significant shifts in their thermostability. HDAC1 and HDAC2 are not detected by the Tm-based approach of (6), because the higher variance between the replicates of the panobinostat-treated condition leads to them being eliminated by the filter heuristics of that method (Table III). In contrast, NPARC naturally takes the variance into account in the computation of the F-statistic and does not require such filtering steps. F–H, Data and curve fits for known non-HDAC targets.
Beyond Two-group Comparisons
Because NPARC is based on analysis of variance (ANOVA), it admits experimental designs in which the covariate has multiple levels. An example is the data set for the BCR-ABL inhibitor dasatinib, which comprises measurements on cells treated at two different concentrations as well as untreated cells. NPARC successfully identified known targets of dasatinib (supplemental Fig. S4).
Replicate Agreement and Model Fit Diagnostics
Application of the Tm -based approach by (6) to the panobinostat data failed to detect HDAC1 and HDAC2. This was because the data for these proteins had relatively high variance in the drug-treated condition, as is visible in Fig. 3A–3B. This led to their exclusion according to one of the data quality filter criteria of that method (Table III), namely the criterion that asks for sufficiently high coefficients of determination (R). In contrast, a better and statistically sound trade-off between variability and effect size is an integral part of NPARC and does not require an ad hoc filter criterion.To further assess the price of the various filter criteria of the Tm-based approach by (6), we tabulated the numbers of proteins affected by them in each of the five data sets. These proteins would, in principle, not be detectable by that method, no matter how strong the effect. Their numbers amounted to 14–25% of the total numbers of proteins for which melting points could be determined in both replicates (Fig. 1E and Table IV) and to 21–32% of all proteins irrespective of melting point availability (supplemental Fig. S2). In contrast, the F-test of NPARC could be applied to all proteins irrespective of these or similar criteria, a fact which contributed to the increased protein coverage and sensitivity of NPARC.
Table IV
Coverage of proteins applicable for hypothesis testing by the original TPP analysis workflow [6]
Data set
Tm outside measured range
Tm available but curves not passing a priori filters
Tm available and curves passing a priori filters
ATP data
220
1004
2953
Dasatinib 0.5 μm data
689
768
3168
Dasatinib 5 μm data
667
507
2980
Panobinostat data
320
461
2868
Staurosporine data
621
631
3253
Effects Beyond Those on the Melting Point
Many of the proteins detected by NPARC displayed reproducible changes in curve shape, whereas their Tm-shifts were small, and not considered significant by the Tm-based approach (Fig. 4). An example is the effect of staurosporine on protein kinase C beta (PRKCB), shown in Fig. 1C. PRKCB is part of the PKC family, whose members were the first reported staurosporine targets (27, 4) and also exhibit similar characteristics (supplemental Fig. S5).
Fig. 4.
NPARC is sensitive to small but reproducible The plot compares the effect size measure used by NPARC, namely RSS0 − RSS1 (y axis), to the Tm-difference (x axis) for those proteins in the staurosporine data set for which Tm estimates could be obtained. Proteins with Benjamini-Hochberg adjusted p values ≤ 0.01 are marked in red if they were exclusively found by NPARC, and in green if they were also detected by the Tm-based approach of (6). NPARC detects targets with small Tm-differences if the measurements are reproducible between replicates.
NPARC is sensitive to small but reproducible The plot compares the effect size measure used by NPARC, namely RSS0 − RSS1 (y axis), to the Tm-difference (x axis) for those proteins in the staurosporine data set for which Tm estimates could be obtained. Proteins with Benjamini-Hochberg adjusted p values ≤ 0.01 are marked in red if they were exclusively found by NPARC, and in green if they were also detected by the Tm-based approach of (6). NPARC detects targets with small Tm-differences if the measurements are reproducible between replicates.Further examples include the effects of staurosporine on RanGTP binding tRNA export receptor exportin-T (XPOT) and two members of the p38 MAPK signaling pathway: Mitogen-activated protein kinase 14 (MAPK14) and MAP kinase-activated protein kinase 2 (MAPKAPK2) (Fig. 4); and the effect of dasatinib on Bruton tyrosine kinase (BTK), an important drug target in B-cell leukemia (Fig. 1B).
Missing Melting Point Estimates
For highly thermostable proteins, the Tm in one or more of the treatment conditions can be outside of the tested temperature range of a TPP experiment (Fig. 1E). One example is NAD(P)H quinone dehydrogenase 2 (NQO2), a cytosolic flavoprotein and a common off-target of kinase inhibitors (28, 29, 30). In concordance with previous CETSA studies that found NQO2 to be highly stable (31), we observed denaturation only beginning at 67 °C (Fig. 1D). Staurosporine treatment further stabilized NQO2 to an extent that it showed no sign of melting in the tested temperature range. The Tm -based approach by (6) will discard such proteins, in order to avoid potential problems from extrapolation of the fit beyond the measured temperature range. In contrast, the functional data analysis approach of NPARC can detect changes in any part of the melting curves, without reference to a single point such as Tm.
Sensitivity and Specificity
So far, we have described increased sensitivity of NPARC, i.e. its ability to detect more true targets. However, this is only useful if at the same time specificity is maintained, i.e. if false positive detection remains under control. To compare these performance characteristics between NPARC and the Tm-based approach, we computed pseudo receiver operator characteristic (ROC) curves for each of these methods on the staurosporine data and the ATP data, using as pseudo ground truth lists of expected targets from Gene Ontology annotation (Table II). Here, the term pseudo refers to the fact that these target lists, and hence the ROC curves, are only approximations of the truth; however, the relative ranking of two methods in such a pseudo-ROC comparison is likely to be faithful even in the presence of such approximation error (32).Fig. 5 shows the results of NPARC on both data sets, as well as those of the z-test of the Tm-based approach by (6) applied to the individual replicates (displayed as continuous lines parameterized by the z-cutoff), and those of the full procedure of (6) (shown by isolated points, because of its single, fixed cutoff). On the staurosporine data, the full procedure of (6) performs close to NPARC. For the individual z-tests, as well as overall on the ATP data, NPARC shows superior performance. Given that the decision rule set of (6) (listed in Table V) and its cutoff parameters were developed and tuned partly on the staurosporine data, these results indicate that NPARC has fewer “fudge parameters” and is likely to be superior in applications to new data sets.
Fig. 5.
Sensitivity and specificity. Shown are pseudo-ROC curves, with expected hits (as a proxy for true positives) along the y axis and unexpected hits (as a proxy for false positives) along the x axis. The curves are obtained by varying the p value cutoff of the F-test of NPARC (which is computed across replicates), and of the z-test of the Tm-based approach (6) (which is computed separately for each replicate). The asterisks indicate the result from the decision rules of (6) on the z-test results (Table V). The dots indicate a threshold of 0.01 on the Benjamini-Hochberg adjusted p values from NPARC (derived on both replicates in parallel) and on the Benjamini-Hochberg adjusted p values from the Tm-based approach (computed individually for each replicate). NPARC is modestly better than the Tm-based approach on the staurosporine data (A), and substantially better on the ATP data (B). The proteins found by NPARC at Benjamini-Hochberg adjusted p value ≤ 0.01 are also shown in supplemental Figs. S6 and S7.
Table V
Decision ruleset applied in the original TPP analysis workflow [6] to combine z-test p values across replicates in an experimental design with two biological replicates
Rule number
Rule
1
The Benjamini-Hochberg adjusted z-test p values fulfill predefined thresholds in each replicate.
2
Both Tm differences are either positive or negative in the two biological replicates.
3
The smallest absolute difference between treatment and vehicle Tm is greater than the absolute Tm difference between the two vehicle experiments.
Sensitivity and specificity. Shown are pseudo-ROC curves, with expected hits (as a proxy for true positives) along the y axis and unexpected hits (as a proxy for false positives) along the x axis. The curves are obtained by varying the p value cutoff of the F-test of NPARC (which is computed across replicates), and of the z-test of the Tm-based approach (6) (which is computed separately for each replicate). The asterisks indicate the result from the decision rules of (6) on the z-test results (Table V). The dots indicate a threshold of 0.01 on the Benjamini-Hochberg adjusted p values from NPARC (derived on both replicates in parallel) and on the Benjamini-Hochberg adjusted p values from the Tm-based approach (computed individually for each replicate). NPARC is modestly better than the Tm-based approach on the staurosporine data (A), and substantially better on the ATP data (B). The proteins found by NPARC at Benjamini-Hochberg adjusted p value ≤ 0.01 are also shown in supplemental Figs. S6 and S7.
DISCUSSION
Thermal proteome profiling offers the possibility to comprehensively characterize ligand-protein interactions on a proteome-wide scale in living cells. However, the method poses the analytical challenge of how to identify statistically significant shifts in thermostability among thousands of measurements.To address this challenge, we introduced a functional data analysis approach to test for treatment effects by comparing competing models by their goodness-of-fit. This enables detection of treatment effects even if a (de-)stabilization of a protein is not captured by a single summary parameter like the Tm. The presented method is based on a sound statistical foundation and does not rely on hard-to-choose cutoff or tuning parameters. We showed that our method compares favorably to previous approaches with respect to sensitivity and specificity for several exemplary data sets, including ones with specific and ones with promiscuous binders.The approach fits into the framework of analysis of variance (ANOVA) or linear models and can thus be extended to experimental designs more complex than treatment-control comparisons, such as multiple levels (e.g. drug concentrations) per covariate, multiple covariates and interactions.The suggested framework is flexible regarding the mean function used to represent the melting behavior and can be adapted to the particular biological process of interest. To represent nonlinear relationships, approaches include locally linear regression (33), spline regression (34, 35) and nonlinear parametric regression. Here, we chose the latter as it incorporates a priori knowledge about the data and thus has favorable estimation efficiency. For example, sigmoid curves have horizontal asymptotes at both sides of the temperature range. In contrast, splines and local regression tend to overfit data near the boundaries of the observation range.In a cellular environment we occasionally observe nonsigmoid melting curves for subsets of proteins. One possible reason is the presence of protein subpopulations each with distinct melting curves (16). For example, the formation of protein complexes, the binding to other molecules, or the localization in cellular compartments can lead to deviations from the idealized sigmoid melting curve expected from the same protein in purified form. Our model currently does not account for such systematic and reproducible shape deviations. This could be adapted in future work by adding a low-parametric systematic modification to the sigmoid mean function.We have considered CETSA experimental designs, where the temperature is the major experimental variable and drug concentration is either zero or a chosen value. It appears relatively straightforward to extend NPARC to the isothermal dose response (ITDR) design (6, 7) where temperature is held constant and the drug concentration is varied across a range of values. A further extension of interest would be to 2D-TPP (11) where both factors are changed.We employ the same “average” null distribution for all proteins, which we obtain by estimating its parameters (d1, d2, σ0) from the distributions of residuals across all proteins. It is conceivable that determining null distributions in a protein dependent manner, for instance by stratification, could increase the overall power of the method.The here presented approach is likely to increase the accuracy of profiling protein-ligand interactions in living cells. We provide an open source R package NPARC, and all computations reported in the figures and tables of this article can be reproduced by running the Rmarkdown script provided in reference (36).
Additional Files
The following Figs. and Tables can be found in the Supplementary MaterialAdditional file 1 — Supplementary MethodsDetailed description of the fitting procedures for the scaling parameter, melting points, and mean functions of the model.Additional file 2 — Table S1Spreadsheet containing the results of the NPARC approach and of the Tm-based approach for all data sets listed in Table I.Additional file 3 — Supplemental Figs. S1–S5PDF containing Supplemental Figs. S1–S5Additional file 4 — Supplemental Fig. S6All proteins detected by the NPARC approach with Benjamini-Hochberg adjusted F-test p values ≤ 0.01 in the staurosporine data.Additional file 5 — Supplemental Fig. S7All proteins detected by the NPARC approach with Benjamini-Hochberg adjusted F-test p values ≤ 0.01 in the ATP data.
DATA AVAILABILITY
Datasets used in this work were downloaded as spreadsheets from journal websites. Table 1 summarizes all datasets used in this work and lists all respective references. The TPP-TR experiments based on staurosporine-, panobinostat-, ATP-, and dasatinib treatments are included in the supplemental materials of references (5, 6, 9). All results generated from this data are provided in the supplemental material attached to this work. A vignette describing the workflow and relevant code to reproduce the results is available for download (36). An open-source R package is available from https://github.com/Huber-group-EMBL/NPARC, and provision via Bioconductor is intended.
Authors: Holger Franken; Toby Mathieson; Dorothee Childs; Gavain M A Sweetman; Thilo Werner; Ina Tögel; Carola Doce; Stephan Gade; Marcus Bantscheff; Gerard Drewes; Friedrich B M Reinhard; Wolfgang Huber; Mikhail M Savitski Journal: Nat Protoc Date: 2015-09-17 Impact factor: 13.491
Authors: Friedrich B M Reinhard; Dirk Eberhard; Thilo Werner; Holger Franken; Dorothee Childs; Carola Doce; Maria Fälth Savitski; Wolfgang Huber; Marcus Bantscheff; Mikhail M Savitski; Gerard Drewes Journal: Nat Methods Date: 2015-11-02 Impact factor: 28.547
Authors: Kenneth M Comess; Shaun M McLoughlin; Jon A Oyer; Paul L Richardson; Henning Stöckmann; Anil Vasudevan; Scott E Warder Journal: J Med Chem Date: 2018-06-06 Impact factor: 7.446
Authors: André Mateus; Jacob Bobonis; Nils Kurzawa; Frank Stein; Dominic Helm; Johannes Hevler; Athanasios Typas; Mikhail M Savitski Journal: Mol Syst Biol Date: 2018-07-06 Impact factor: 11.429
Authors: Kilian V M Huber; Karin M Olek; André C Müller; Chris Soon Heng Tan; Keiryn L Bennett; Jacques Colinge; Giulio Superti-Furga Journal: Nat Methods Date: 2015-09-21 Impact factor: 28.547
Authors: Mikhail M Savitski; Nico Zinn; Maria Faelth-Savitski; Daniel Poeckel; Stephan Gade; Isabelle Becher; Marcel Muelbaier; Anne J Wagner; Katrin Strohmer; Thilo Werner; Stephanie Melchert; Massimo Petretich; Anna Rutkowska; Johanna Vappiani; Holger Franken; Michael Steidel; Gavain M Sweetman; Omer Gilan; Enid Y N Lam; Mark A Dawson; Rab K Prinjha; Paola Grandi; Giovanna Bergamini; Marcus Bantscheff Journal: Cell Date: 2018-03-15 Impact factor: 41.582
Authors: Alice L Herneisen; Saima M Sidik; Benedikt M Markus; David H Drewry; William J Zuercher; Sebastian Lourido Journal: ACS Chem Biol Date: 2020-07-06 Impact factor: 5.100
Authors: Fraser D Johnson; John Ferrarone; Alvin Liu; Christina Brandstädter; Ravi Munuganti; Dylan A Farnsworth; Daniel Lu; Jennifer Luu; Tianna Sihota; Sophie Jansen; Amy Nagelberg; Rocky Shi; Giovanni C Forcina; Xu Zhang; Grace S W Cheng; Sandra E Spencer Miko; Georgia de Rappard-Yuswack; Poul H Sorensen; Scott J Dixon; Udayan Guha; Katja Becker; Hakim Djaballah; Romel Somwar; Harold Varmus; Gregg B Morin; William W Lockwood Journal: Cell Rep Date: 2022-02-08 Impact factor: 9.423
Authors: Simone Di Sanzo; Katrin Spengler; Anja Leheis; Joanna M Kirkpatrick; Theresa L Rändler; Tim Baldensperger; Therese Dau; Christian Henning; Luca Parca; Christian Marx; Zhao-Qi Wang; Marcus A Glomb; Alessandro Ori; Regine Heller Journal: Nat Commun Date: 2021-11-18 Impact factor: 14.919
Authors: Nils Kurzawa; Isabelle Becher; Sindhuja Sridharan; Holger Franken; André Mateus; Simon Anders; Marcus Bantscheff; Wolfgang Huber; Mikhail M Savitski Journal: Nat Commun Date: 2020-11-13 Impact factor: 14.919
Authors: Marcin Luzarowski; Rubén Vicente; Andrei Kiselev; Mateusz Wagner; Dennis Schlossarek; Alexander Erban; Leonardo Perez de Souza; Dorothee Childs; Izabela Wojciechowska; Urszula Luzarowska; Michał Górka; Ewelina M Sokołowska; Monika Kosmacz; Juan C Moreno; Aleksandra Brzezińska; Bhavana Vegesna; Joachim Kopka; Alisdair R Fernie; Lothar Willmitzer; Jennifer C Ewald; Aleksandra Skirycz Journal: Commun Biol Date: 2021-02-10