← Bioinformatics Series

Clinical Evidence

A waterfall plot asks who responded, a survival curve asks who lived longer, and a forest plot asks whether an effect holds up from one study to the next.

The earlier parts of this series followed one question through molecular data: how a tumor ends up with p21, the protein made from CDKN1A, switched off. Part 1 laid out three routes (a TP53 mutation, an MDM2 amplification, a methylated CDKN1A promoter), and Part 6 recovered them as clusters of patients and left one question for this part: do the clusters differ in outcome? The 100 patients in Figures 1 to 4 are the same simulated patients as in those two parts, so their routes and clusters carry over. Figures 5 and 6 pool eight separate simulated studies.

Waterfall plots

In oncology, a waterfall plot usually means one bar per patient, sorted by the best percent change in target-lesion size during treatment. Bars below zero are tumors that shrank and bars above are tumors that grew. The plot keeps every patient visible where a response rate would collapse them into one number.

Rank by
Color by

Figure 1. The 100 simulated patients as ranked bars. Ranked by best response, each bar is the best percent change from baseline in the summed diameters of the target lesions. The −30% line is the partial-response threshold. The +20% line marks progression, which RECIST measures from the smallest sum on study (the nadir); for a best change the nadir is the baseline. A dot marks a new lesion, which counts as progression whatever the bar does. A complete response means every target lesion has gone (−100%). Tab to the chart and use the arrow keys, or hover or tap a bar, to read one patient.

Colored by TP53 status, the mutant tumors sit toward the growth end: their median best change is , against for wildtype tumors. Colored by Part 6's clusters, two more things show. The MDM2-amplified tumors are TP53-wildtype but respond almost as poorly as the mutants (median ), because MDM2 keeps p53 off without touching the gene. The methylated tumors respond best of all (median ), although their CDKN1A is as low as the mutants'. In this simulation, response was made to depend on whether p53 itself works, not on p21, so a TP53 test alone would misplace both groups.

Ranked bars work for molecular quantities too, as long as the plot names what it ranks. Rank by CDKN1A and the TP53-mutant bars collect at the low end, with the MDM2 and methylated tumors among them: three routes to the same low p21. Rank by tumor mutation burden (TMB, mutations per megabase of sequenced DNA) and the mutants move to the high end, with a median of against . Real cohorts show the same association, and it has more than one explanation. Without p53, cells with damaged DNA keep dividing where they would otherwise stop or die, so mutations pile up. Mutagens such as tobacco smoke also raise TMB and cause TP53 mutations at the same time. In this cohort the association was built in.

Survival analysis

A waterfall looks at one moment. Survival analysis follows patients over time. The Kaplan-Meier curve estimates, at each time, the fraction of a group who have not yet had the event, here death:

\[ S(t) = \prod_{i: t_i \le t} \left( 1 - \frac{d_i}{n_i} \right) \]

where \( d_i \) is the number of deaths at time \( t_i \) and \( n_i \) is the number of patients still at risk just before it. Each step down is a death. Tick marks on the curve are censored patients, who were lost to follow-up or still alive when the study ended. They are known to have survived at least that long.

Figure 2. One Kaplan-Meier curve for the 100 simulated patients: deaths and censored, median survival months. The table counts patients still at risk. The shaded band is the 95% confidence interval (Greenwood's variance on the log(−log) scale). It tracks how many patients remain: wide at 12 months with at risk, at 48 months with .

Without censoring, patients still alive at the end would have to be dropped, which wastes what is known about them and can bias the curve, or counted as deaths, which makes survival look worse than it is. The Kaplan-Meier estimator keeps each patient in the risk set until their last follow-up and then removes them from later denominators. That is only fair if censoring is non-informative: the patients who leave follow-up must have the same future risk as those who stay. Here censoring comes from staggered enrolment and random loss to follow-up, both independent of each patient's risk.

Splitting patients by a marker asks whether it goes with survival. Figure 3 splits by CDKN1A expression at the median, a cut fixed before looking. GAPDH, a housekeeping gene that the simulation gave no effect on anything, is the control. The third option is Part 6's four clusters.

Split by

Figure 3. Kaplan-Meier curves for the same 100 simulated patients, split by expression of one gene (red: at or above the cut, blue: below) or by Part 6's consensus clusters, each named after the planted route most of its members took. The p-value is valid for one cut chosen in advance. Figure 4 shows what moving the slider in search of a better cut does to it.

A gap between curves is not yet a difference. Curves can cross, and a tail with a handful of patients at risk can swing by a lot. The log-rank test compares observed with expected deaths in each group at every death time, and asks how surprising the total gap would be if the groups shared one survival curve. A small p says the curves differ, not by how much. The hazard ratio from a Cox model gives the size: how many times faster one group dies than the other at any moment.

Split at the median, low CDKN1A goes with shorter survival (p = , hazard ratio ). GAPDH does not (p = ). The four clusters do not separate at 0.05 (p = ), although three of them are routes to low p21. A four-group test spreads its evidence over three comparisons, and the MDM2 cluster has only patients. Pooling the three routes against the intact cluster gives p = , still short of 0.05, partly because of the intact cluster's tumors took a route in the simulation; clustering placed them by their data, which looked intact enough. The simulation's truth is simpler than either view: survival depends on CDKN1A level and stage, and the routes matter only through how far they push CDKN1A down.

Choosing the cut

The median is not the only possible cut, and it is tempting to try others until the curves separate best. Altman, Lausen, Sauerbrei and Schumacher (1994) showed what that does to the p-value. The smallest p over cuts between the 10th and 90th percentiles is far smaller than a single test would give, and a marker with no effect at all reaches p < 0.05 much more often than one time in twenty. Figure 4 runs that experiment on this cohort. It does the log-rank test at every cut from the 10th to the 90th percentile, then shuffles the gene's values across patients times, which breaks any link to survival, and repeats the whole scan each time.

Scan

Figure 4. Cutpoint scan for the 100 simulated patients. The colored line is −log10 of the log-rank p at each cut. The gray lines are 30 of the scans after shuffling the gene's values. The dashed line is p = 0.05; the solid line is the 5% point of the shuffled scans' smallest p, the bar a scanned result has to clear. The marker is Figure 3's cut: tab to the chart and use the arrow keys, or click a percentile, to move it.

GAPDH dips below 0.05 at of the 81 cuts, best at the percentile (p = ). A paper that reported only that cut would present a prognostic housekeeping gene. Against the shuffled scans the same minimum is unremarkable: of them reach one at least as small. CDKN1A clears the corrected bar, with a corrected p of . The 5% point of the shuffled minima, about , lands where Altman and colleagues' correction formula puts it: a smallest p of 0.002 over this range corresponds to about . The safe options are to fix the cut before looking, to correct the minimum p for the scan, or to test the chosen cut in a second cohort.

The same arithmetic applies across genes. Test 20,000 genes for survival and about 1,000 will reach p < 0.05 by chance, which is why survival screens need the same multiple-testing correction as a differential expression analysis.

Pooling across studies

A single cohort can mislead, so a meta-analysis pools results across independent studies to ask whether an effect holds in other tissues and labs. In a forest plot each row is a study, drawn as an effect size with its confidence interval. The square marks the estimate and its area is proportional to the study's weight. The diamond at the bottom is the pooled estimate.

Pooling model

Figure 5. Eight simulated studies, separate from the cohort above, each reporting the log2 fold change of CDKN1A in TP53-mutant against wildtype tumors. The dashed line is no change. Click a row, or tab to it and press Enter, to leave that study out. Leaving out , the heaviest study under fixed effect, moves the random-effects estimate from to . Leaving out does more to the spread: I² falls from to .

The two pooling models answer different questions. Fixed effect assumes every study measures one true effect and weights each by the inverse of its variance, so the large studies dominate: weights run from to . Random effects assumes the true effect varies from study to study. It estimates that between-study variance, τ² (by DerSimonian and Laird's method, here ), and adds it to every study's own variance. Because τ² is the same for all studies it flattens the weights, to to , and it widens the diamond.

The I² statistic estimates what share of the spread between studies comes from real differences rather than sampling error. Here it is , which calls for an explanation of why the studies differ: tissue, assay platform, who was enrolled, or how the data were analyzed. Leaving studies out one at a time, as the rows of Figure 5 allow, shows whether a single study is carrying the conclusion. Here none is: every leave-one-out interval stays below zero ().

Journals tend to prefer significant results, which can bias the published record. The funnel plot draws each study's effect against its standard error, with the most precise studies at the top. With one true effect and no bias, studies scatter symmetrically inside a cone that narrows toward the top. With only a few studies, asymmetry is a reason to look closer and does not establish bias.

Figure 6. Funnel plot of the studies in Figure 5. It follows that figure's pooling model and exclusions; left-out studies are hollow. The vertical line is the pooled estimate, and the cone is that estimate ± 1.96 standard errors, widened by τ² under random effects. Here the asymmetry runs against the publication-bias story: the two smallest studies, , are the only two whose intervals cross zero, and the simulation gave those tissues weaker true effects. Egger's regression intercept is ; bias toward publishing strong negative effects would push it below zero.

A funnel can be lopsided for reasons that have nothing to do with which studies got published: small studies may enroll different patients, use different assays, or, as here, come from tissues where the effect really is weaker. Sterne and colleagues (2011) list heterogeneity among the first causes to rule out before reading asymmetry as bias.

Back to p21

Part 1 laid out the routes to a silent CDKN1A, Part 4 asked whether the R175H mutant still binds the CDKN1A promoter, Part 5 looked for CDKN1A among the genes that fall when p53 is broken, and Part 6 recovered the routes as clusters. Here, low CDKN1A goes with shorter survival at a cut fixed in advance (p = ) and still after correcting for a scan of every cut (p = ).

What these figures cannot say is why. TP53 status on its own gives p = , weaker than CDKN1A, though TP53 sits at the start of the chain. In a Cox model that also includes CDKN1A, TP53's hazard ratio is (p = ), and adding Part 6's clusters to a model with CDKN1A does not reach significance (likelihood-ratio p = ), as it should not: in the simulation the clusters have no effect of their own. That pattern fits CDKN1A as the cause, but it fits CDKN1A as a marker for something else equally well. Only the simulation knows that survival was written as a function of CDKN1A and stage. In real data, a survival split shows that a marker and an outcome travel together, and a meta-analysis shows whether that holds across studies. Neither shows that changing the marker would change the outcome.