## Graphical Abstract

**Highlights**

Method for the analysis of response curves from thermal proteome profiling (TPP).

NPARC uses nonparametric statistics and provides false discovery-rate (FDR) control.

Increased proteome coverage and sensitivity to identify drug-binding proteins.

## 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.

- Drug targets
- algorithms
- biostatistics
- tandem mass spectrometry
- mathematical modeling
- functional data analysis
- 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)^{1} 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 (*T*_{m}), which is defined as the temperature of half-maximum relative abundance (horizontal line in Fig. 1*A*). 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 *T*_{m} (Fig. 1*B*–1*C*). Second, the true *T*_{m} of a protein can lie outside of the measured temperature range, which impairs its estimation (Fig. 1*D*). Both scenarios can result in important targets being missed in the analysis (Fig. 1*E*). The third reason is a more subtle statistical one: hypothesis tests using only the point estimates of *T*_{m} 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.

Here, we propose an alternative approach that compares whole curves instead of summary parameters and does not rely on *T*_{m} 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 *T*_{m}-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 *T*_{m}-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.

##### 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. 2*A*). The alternative model posits two condition-specific functions: μ* _{T}*(

*t*) for the treatment condition and μ

*(*

_{V}*t*) for the vehicle condition (Fig. 2

*B*). 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

*x*is the measured value at temperature

_{t,i,c}*t*for experimental replicate

*i*and condition

*c*∈ {

*V*,

*T*}, and the summations extend over all temperatures, replicates, and conditions.

##### Choice of the Mean Function

The mean functions μ_{0}(*t*), μ* _{T}*(

*t*) and μ

*(*

_{V}*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 prescription

The 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
*d*_{2}/*d*_{1} > 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 μ* _{0}(t*). If the residuals were independent and identically normal distributed, this distribution would be given by an analytical formula, namely that of the

*F*(

*d*

_{1},

*d*

_{2})-distribution with parameters

*d*

_{1},

*d*

_{2}> 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

*RSS*

_{0}and

*RSS*

_{1}. 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*(

*d*

_{1},

*d*

_{2})-distributions is quite flexible, and we can approximate the distribution of the

*F*-statistic (Eq. 4) on data occurring in practice with an

*F*(

*d*

_{1},

*d*

_{2})-distribution with different “effective degrees of freedom”

*d*

_{1},

*d*

_{2}. To this end, we separately approximate the numerator and denominator of

*F*as

*-distributed random variables χ*

^{2}

^{2}(d_{1}), χ

^{2}(d_{2}) has an

*F*(

*d*

_{1},

*d*

_{2})-distribution (24). The scale parameter σ

_{0}

^{2}and the effective degrees of freedom

*d*

_{1}and

*d*

_{2}are estimated from the empirical distributions—across proteins—of

*RSS*

_{0}and

*RSS*

_{1}. Thus, we assume that σ

_{0}

^{2},

*d*

_{1},

*d*

_{2}are the same for all proteins. We estimate σ

_{0}

^{2}from the moments of

*RSS*

_{0}−

*RSS*

_{1}as

*RSS*

_{0}−

*RSS*

_{1}(see Supplementary Methods for details). Then,

*d*

_{1}and

*d*

_{2}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 *d*_{1}, *d*_{2} 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. 2*C*.

## 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. 3*A*–3*E*) 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. 3*F*–3*H*). 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. 1*B*). 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.

##### 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 *T*_{m} -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. 3*A*–3*B*. 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 ^{2}*). 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 *T*_{m}-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. 1*E* 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.

##### Effects Beyond Those on the Melting Point

Many of the proteins detected by NPARC displayed reproducible changes in curve shape, whereas their *T*_{m}-shifts were small, and not considered significant by the *T*_{m}-based approach (Fig. 4). An example is the effect of staurosporine on protein kinase C beta (PRKCB), shown in Fig. 1*C*. 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).

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. 1*B*).

##### Missing Melting Point Estimates

For highly thermostable proteins, the *T*_{m} in one or more of the treatment conditions can be outside of the tested temperature range of a TPP experiment (Fig. 1*E*). 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. 1*D*). Staurosporine treatment further stabilized NQO2 to an extent that it showed no sign of melting in the tested temperature range. The *T*_{m} -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 *T*_{m}.

##### 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 *T*_{m}-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 *T*_{m}-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.

## 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 *T*_{m}. 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 (*d*_{1}, *d*_{2}, σ_{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 Material

Additional file 1 — Supplementary Methods

Detailed description of the fitting procedures for the scaling parameter, melting points, and mean functions of the model.

Additional file 2 — Table S1

Spreadsheet containing the results of the NPARC approach and of the *T*_{m}-based approach for all data sets listed in Table I.

Additional file 3 — Supplemental Figs. S1–S5

PDF containing Supplemental Figs. S1–S5

Additional file 4 — Supplemental Fig. S6

All 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. S7

All 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.

## Acknowledgments

We thank Sindhuja Sridharan (EMBL Heidelberg) for help with figure preparation. We thank EMBL and GlaxoSmithKline for supporting the work.

## Footnotes

Author contributions: D.C., K.B., H.F., S.A., M.B., M.M.S., and W.H. designed research; D.C., K.B., H.F., and N.K. performed research; D.C., K.B., H.F., M.M.S., and W.H. analyzed data; D.C., K.B., H.F., S.A., N.K., M.B., M.M.S., and W.H. wrote the paper.

↵* K.B. is supported by a Cambridge Cancer Centre studentship. S.A. is funded by the Deutsche Forschungsgemeinschaft, SFB 1036. W.H. acknowledges funding from the European Commission's H2020 Programme, Collaborative research project SOUND (Grant Agreement no 633974). One or more authors has an actual or perceived conflict of interest with the contents of this article. H.F., M.S., D.C., and M.B. are employees and/or shareholders of GlaxoSmithKline. N.K. is supported by a fellowship of the European Molecular Biology Laboratory International PhD Program.

↵

^{}This article contains supplemental material.↵

^{1}The abbreviations used are:- TPP
- Thermal proteome profiling
- CETSA
- Cellular thermal shift assay
- FDR
- False discovery rate
- H
_{0} - Null hypothesis
- NPARC
- Nonparametric analysis of response curves
- ROC
- Receiver operating characteristic
- RSS
- Residual sum of squares
*T*_{m}- Melting point.

- Received April 8, 2019.
- Revision received September 1, 2019.

- © 2019 Childs et al.

*Author's Choice*—Final version open access under the terms of the Creative Commons CC-BY license.