From Molecular Profiles to Drug Response
A pharmacologic response analysis asks whether molecular measurements obtained before treatment are associated with how strongly a cell line or patient responds to a drug. The molecular feature might be gene expression, protein abundance, or a metabolite level. The response might come from a laboratory dose-response assay or from a clinical measurement. The aim is often to identify candidate response markers: features associated with variation in response that could merit independent investigation.
The analysis is not automatically a test of whether changing a molecular feature would change drug sensitivity. An association can arise because the feature tracks another biological property, treatment dose, or assay batch. The experimental unit also matters: multiple wells from one cell line are not necessarily independent biological observations. These considerations should shape the model and the strength of the conclusion.
Choose and Interpret the Response Measurement
A dose-response experiment measures an outcome at several concentrations. For example, the outcome may be the fraction of viable cells, where larger values indicate less drug effect. A summary such as the normalized area under the viability curve (AUC) compresses that curve into one number. With this convention, a larger AUC means greater viability across the tested range and therefore lower sensitivity to the drug. Other response metrics have different directions: a larger maximum inhibition can indicate greater sensitivity, while a larger IC50 generally indicates that more drug is needed for a given inhibitory effect. Always state the metric and its direction.
AUC depends on the concentration range and how the curve is summarized. It is not a universal measure of response independent of assay design. If each biological unit has a single outcome at one administered dose, dose may need to be modeled as a covariate. If response is measured at multiple doses and summarized as AUC, those doses define the response curve; they should not be treated as separate independent biological units in a marker regression.
Worked Example: Summarizing a Viability Curve
Suppose a hypothetical cell line has normalized viability values 0.90, 0.60, and 0.20 at three equally spaced log-concentration values, 0, 1, and 2. Use the trapezoidal rule and divide by the concentration-range width to obtain normalized AUC.
The two trapezoid areas are \(1(0.90+0.60)/2=0.75\) and \(1(0.60+0.20)/2=0.40\). The total area is \(0.75+0.40=1.15\), over a log-concentration range of \(2-0=2\).
This is a dimensionless normalized AUC on the stated viability scale. Under the convention used here, larger values represent greater viability over the tested range, so 0.575 indicates more resistance than a lower AUC would, all else equal. It does not by itself establish a clinically relevant dose or response.
Model an Adjusted Association
For a continuous response and a continuous baseline feature, a basic model can adjust for dose and assay batch. Let \(D_i\) denote dose or exposure when it is relevant to the response definition, and let \(B_i\) denote batch. For a categorical batch with a reference category, a linear model for feature \(k\) is
Here \(\beta_k\) is the adjusted change in mean response associated with a one-unit increase in feature \(k\), comparing observations with the same modeled dose and batch. The coefficient's units are response units per feature unit. If the feature is standardized, its unit is one sample standard deviation; that choice should be reported. If response is normalized AUC, \(\beta_k\) is in normalized AUC units per feature unit.
The usual linear-model inference assumes a correctly specified conditional mean, independent biological-unit errors, and constant error variance for the standard ordinary least-squares standard errors. Exact small-sample \(t\)-tests additionally assume normally distributed errors conditional on the predictors. The predictors must have enough independent variation to estimate separate coefficients: if feature values are completely determined by dose or batch, their adjusted effects cannot be disentangled. Independence concerns biological units, not merely rows in a file.
Worked Example: Estimating a Feature Association With Batch Adjustment
Consider hypothetical normalized AUC values from eight independent cell lines. A baseline feature \(Z\) is recorded in arbitrary units, and each batch contains four cell lines. The table lists the response values in each batch.
| Batch | Feature \(Z\) | Normalized AUC \(Y\) |
|---|---|---|
| A | 0 | 0.28 |
| A | 2 | 0.52 |
| A | 0 | 0.32 |
| A | 2 | 0.48 |
| B | 0 | 0.28 |
| B | 2 | 0.52 |
| B | 0 | 0.32 |
| B | 2 | 0.48 |
Fit \(Y_i=\beta_0+\beta_k Z_i+\delta B_i+\varepsilon_i\), using batch A as the reference. Each batch has mean feature value 1 and mean response 0.40. The within-batch fitted line has intercept 0.30 and slope 0.10: at \(Z=0\), the mean response is 0.30; at \(Z=2\), it is 0.50. Thus the fitted value for each batch is 0.30, 0.50, 0.30, 0.50 in the row order shown. The batch coefficient is 0 in these data.
The residuals in each batch are \(-0.02, 0.02, 0.02, -0.02\). Therefore,
There are three fitted coefficients: an intercept, a feature slope, and a batch effect. Hence \(df=8-3=5\). The feature values, centered within batch, have squared deviations summing to 8, so
Under the model assumptions, a two-sided \(t\)-test with 5 degrees of freedom gives \(p\) about 0.0001. A 95% confidence interval using \(t_{0.975,5}\approx2.571\) is \(0.100\pm2.571(0.008944)\), or approximately (0.0770, 0.123). In this hypothetical analysis, higher feature values are associated with higher normalized AUC, and thus lower sensitivity under the stated convention, conditional on batch. The small sample and model assumptions warrant caution despite the small p-value.
Why Adjustment and Design Checks Matter
Including a covariate in a regression is useful only when the data contain information to separate its association from that of the molecular feature. If all high feature values occur in one batch and all low values in another, batch and feature are confounded in the observed design. A regression may then be unable to estimate their effects separately, or may produce highly unstable estimates when the overlap is merely limited. A coefficient labelled “adjusted” does not guarantee that confounding has been solved.
A short justification for the adjusted slope is obtained by removing the effects of the other predictors first. In a linear model, regress the feature on the intercept and adjustment variables, and retain its residuals \(r_i\). Regress the response on those same adjustment variables, and retain residuals \(q_i\). The least-squares coefficient for the feature is the slope of \(q_i\) on \(r_i\):
This follows by expressing the model's feature contribution using the component of the feature orthogonal to the adjustment-variable columns. The numerator measures how response variation left after adjustment covaries with feature variation left after adjustment; the denominator is the remaining feature variation. If every \(r_i=0\), the denominator is zero and the adjusted feature coefficient is not identifiable. If the residualized feature values are very small, the estimate may be highly uncertain.
Worked Example: Recognizing a Non-Identifiable Marker Effect
Suppose a hypothetical study measures eight cell lines, with feature \(Z=0\) in every line from batch A and \(Z=1\) in every line from batch B. The proposed model includes an intercept, \(Z\), and a batch indicator \(B\), where \(B=0\) for A and \(B=1\) for B.
In this design \(Z_i=B_i\) for every observation. Thus \(\beta_k Z_i+\delta B_i=(\beta_k+\delta)B_i\). The data can estimate the combined batch-feature contrast, but cannot determine how much belongs to \(\beta_k\) and how much to \(\delta\). Equivalently, after regressing \(Z\) on batch, every residual \(r_i=0\), so \(\sum_i r_i^2=0\) in the adjusted-slope formula.
The correct conclusion is that the adjusted feature effect is not identifiable from these data. Reporting an ordinary feature-only association as a batch-adjusted marker result would overstate what the design supports. More overlap of feature values across batches, or new data with such overlap, is needed to separate the effects.
Multiple Candidate Features and Practical Diagnostics
Testing many molecular features creates a multiple-testing problem: even if none is associated with response, some small p-values can arise by chance. When controlling the false discovery rate, the Benjamini–Hochberg (BH) procedure orders the \(m\) p-values from smallest to largest, \(p_{(1)}\le\cdots\le p_{(m)}\), and finds the largest rank \(j\) satisfying \(p_{(j)}\le j\alpha/m\). It selects the tests ranked 1 through \(j\); if no rank qualifies, it selects none. This procedure controls the expected proportion of false discoveries among selected tests under independence and under certain forms of positive dependence. The dependence structure of molecular tests should be considered when choosing and describing a procedure.
Worked Example: Applying Benjamini–Hochberg to Three Tests
Suppose three hypothetical marker tests have ordered p-values 0.010, 0.040, and 0.200. At false-discovery-rate level \(\alpha=0.05\), the BH thresholds are \(1(0.05)/3=0.0167\), \(2(0.05)/3=0.0333\), and \(3(0.05)/3=0.0500\).
The first p-value qualifies because \(0.010\le0.0167\). The second does not, because \(0.040>0.0333\), and the third does not, because \(0.200>0.0500\). The largest qualifying rank is 1, so only the first test is selected. The second test is not selected merely because it is the second-smallest p-value.
Diagnostics should examine both the fitted model and the study design. Useful checks include plotting response against the feature within each batch, verifying that feature values overlap across batches and doses, and inspecting residuals for curvature, changing spread, or influential observations. For a dose-response summary, inspect the underlying curves as well: an AUC can conceal unusual curve shapes or a response that differs only outside the measured range. Check whether missing profiles or assay failures systematically exclude particular response groups.
Technical replicate wells can improve measurement precision, but they do not automatically increase the number of independent cell lines or patients. If wells from the same biological unit are analyzed as independent rows, standard errors can be too small because outcomes from those wells may be correlated. Summarizing technical replicates at the biological-unit level, or using a model that represents their dependence, is more appropriate. For patient cohorts, relevant clinical covariates and the way dose or exposure was assigned also require consideration.
Choosing an Approach and Stating Its Limits
| Approach | Useful when | Main caution |
|---|---|---|
| Single-feature regression with covariate adjustment | The response is quantitative and the association of a prespecified feature is the target. | Linear mean structure and adequate feature overlap are needed for straightforward interpretation. |
| Response-curve summary, such as normalized AUC | Each biological unit has measurements across a dose range and one interpretable summary is desired. | The summary depends on the measured range and can conceal curve-shape differences. |
| Separate unadjusted comparisons by feature | Exploratory visualization or a simple description is useful. | Confounding by batch, dose, or other variables is not controlled by the comparison. |
| Flexible or nonlinear response model | Plots or subject-matter knowledge suggest curvature or a dose-by-feature interaction. | More parameters require adequate sample size and careful validation. |
No single model resolves every limitation. A linear feature term may miss nonlinear associations; an additive model may miss a feature whose association changes with dose. Adjustment can increase uncertainty, and it cannot correct a complete lack of overlap or unmeasured confounding. A low p-value is not a measure of predictive usefulness, clinical benefit, or reproducibility. Candidate markers need evaluation in independent data and in a design appropriate to the intended use.
Common Mistakes
- Reversing the response direction. With viability AUC, a higher value means higher viability and usually lower sensitivity; with an inhibition metric, the direction may be the opposite. Define it before interpreting a coefficient.
- Calling an adjusted association causal. Adjustment for measured batch and dose does not guarantee exchangeability or remove unmeasured confounding.
- Treating wells as independent subjects. Replicate wells from one cell line do not provide the same independent information as additional cell lines.
- Assuming adjustment fixes confounding. If feature and batch are perfectly aligned, their separate effects cannot be estimated from that design.
- Applying BH as a rank-by-rank list of rejections. Find the largest rank that meets its threshold; if that rank is \(j\), select ranks 1 through \(j\). If rank 2 fails but rank 1 qualifies, only rank 1 is selected.
Check Your Understanding
Answer each question, then compare your reasoning with the solutions.
- If higher normalized viability AUC indicates greater viability, what does a positive feature slope mean in an adjusted model?
- A feature is exactly equal to the batch indicator for every sample. Can its batch-adjusted coefficient be estimated separately from the batch effect?
- Why can treating repeated wells from each cell line as independent observations make inference misleading?
- For ordered p-values 0.010, 0.040, and 0.200 at \(\alpha=0.05\), how many tests does BH select?
- Does a statistically significant adjusted marker association show that changing the marker would change drug response?
Solutions. (1) Higher feature values are associated with higher normalized AUC, hence greater viability and lower sensitivity under that convention, conditional on modeled covariates. (2) No. The feature is perfectly explained by batch, so its separate adjusted coefficient is not identifiable. (3) Wells from the same biological unit may have dependent outcomes; treating them as independent can understate uncertainty. (4) The thresholds are 0.0167, 0.0333, and 0.0500. Only the first qualifies: \(0.010\le0.0167\), while \(0.040>0.0333\) and \(0.200>0.0500\). BH selects only the first test. (5) No. The result is an adjusted association; causal interpretation requires an appropriate design and additional assumptions.