library(gDRcore)
library(gDRtestData)
library(gDRutils)
library(SummarizedExperiment)
library(BumpyMatrix)
library(data.table)

Overview

The standard gDR pipeline fits dose-response curves using a fixed 4-parameter log-logistic model (fit_SE()). Sometimes you need something different: an alternative curve model, custom synergy metrics, Bayesian estimates, or bespoke pharmacology metrics for your assay.

apply_fit() and apply_fits() let you plug any R function into the gDR SE/MAE pipeline without touching the pipeline internals.

These functions, the fit profile registry and fit_drug_response_metrics() now live in gDRutils, next to the curve-fitting maths they are built on. gDRcore re-exports them for one release cycle, so the code below keeps working as written; call them as gDRutils::apply_fit() in new code.

Key properties:

Works with single-agent, combination, and time-course data.
Each fit function writes to its own named assay — no risk of overwriting native gDR assays ("Metrics", "scores", "excess", …).
Results are idempotent: calling twice with the same fit_source overwrites rather than duplicates.
apply_fits() applies N fit functions in a single BumpyMatrix traversal — efficient when you have several metrics to compute on the same data.

Data setup

We use a small synthetic single-agent SE from gDRtestData.

mae <- gDRutils::get_synthetic_data("finalMAE_small.qs2")
sa_se <- mae[["single-agent"]]
sa_se
#> class: SummarizedExperiment 
#> dim: 10 10 
#> metadata(5): identifiers experiment_metadata Keys fit_parameters
#>   .internal
#> assays(5): RawTreated Controls Normalized Averaged Metrics
#> rownames(10): G00002_drug_002_moa_A_72 G00003_drug_003_moa_A_72 ...
#>   G00010_drug_010_moa_A_72 G00011_drug_011_moa_B_72
#> rowData names(4): Gnumber DrugName drug_moa Duration
#> colnames(10): CL00011_cellline_BA_tissue_x_26
#>   CL00012_cellline_CA_tissue_x_30 ... CL00019_cellline_JB_tissue_z_58
#>   CL00020_cellline_KB_tissue_z_62
#> colData names(4): clid CellLineName Tissue ReferenceDivisionTime

The Averaged assay is the input to custom fit functions. Each cell of the BumpyMatrix contains one data.table per (drug × cell line) pair:

avg_cell <- BumpyMatrix::unsplitAsDataFrame(
  assay(sa_se, "Averaged"),
  row.field = "row", column.field = "column"
)
head(avg_cell[avg_cell$row == avg_cell$row[1] &
              avg_cell$column == avg_cell$column[1], ])
#> DataFrame with 6 rows and 6 columns
#>                      row                 column normalization_type
#>              <character>            <character>           <factor>
#> 1 G00002_drug_002_moa_.. CL00011_cellline_BA_..                 RV
#> 2 G00002_drug_002_moa_.. CL00011_cellline_BA_..                 GR
#> 3 G00002_drug_002_moa_.. CL00011_cellline_BA_..                 RV
#> 4 G00002_drug_002_moa_.. CL00011_cellline_BA_..                 GR
#> 5 G00002_drug_002_moa_.. CL00011_cellline_BA_..                 RV
#> 6 G00002_drug_002_moa_.. CL00011_cellline_BA_..                 GR
#>   Concentration         x      x_std
#>       <numeric> <numeric>  <numeric>
#> 1    0.00100000  0.924967 0.01066130
#> 2    0.00100000  0.944433 0.00805874
#> 3    0.00316228  0.739100 0.01484621
#> 4    0.00316228  0.793100 0.01302152
#> 5    0.01000000  0.435933 0.03071503
#> 6    0.01000000  0.481333 0.03811841

Each cell has four columns:

Column	Description
`normalization_type`	`"GR"` or `"RV"`
`Concentration`	drug concentration (µM)
`x`	averaged normalized response (GR value or relative viability)
`x_std`	standard deviation across replicates

The fit_fn contract

A fit function must satisfy this interface:

fit_fn(avg_dt) -> named list (or single-row data.frame / data.table)

avg_dt is a data.table containing the rows for one combination of slicing column values (by default: one normalization_type) for one (drug × cell line) cell.

The returned named list becomes one row in the output assay. Column names are the names of the list elements. fit_source is stamped automatically by the generic layer — do not include it in the return value.

Minimal example

# The simplest possible fit_fn: compute mean and SD of the response
summary_fn <- function(avg_dt) {
  list(
    x_mean = mean(avg_dt$x, na.rm = TRUE),
    x_sd = sd(avg_dt$x, na.rm = TRUE),
    n = NROW(avg_dt)
  )
}

Single-agent: apply_fit()

Apply the summary function to every (drug × cell line × normalization_type) triplet and write results to a custom assay named "custom_summary".

sa_out <- apply_fit(
  sa_se,
  fit_fn = summary_fn,
  data_type = "single-agent",
  output_assay = "custom_summary",
  fit_source = "demo"
)
assayNames(sa_out)
#> [1] "RawTreated"     "Controls"       "Normalized"     "Averaged"      
#> [5] "Metrics"        "custom_summary"

summary_df <- BumpyMatrix::unsplitAsDataFrame(
  assay(sa_out, "custom_summary"),
  row.field = "row", column.field = "column"
)
head(summary_df)
#> DataFrame with 6 rows and 7 columns
#>                      row                 column    x_mean      x_sd         n
#>              <character>            <character> <numeric> <numeric> <integer>
#> 1 G00002_drug_002_moa_.. CL00011_cellline_BA_..  0.491452  0.221336         9
#> 2 G00002_drug_002_moa_.. CL00011_cellline_BA_..  0.466411  0.214770         9
#> 3 G00003_drug_003_moa_.. CL00011_cellline_BA_..  0.417704  0.273562         9
#> 4 G00003_drug_003_moa_.. CL00011_cellline_BA_..  0.419170  0.256418         9
#> 5 G00004_drug_004_moa_.. CL00011_cellline_BA_..  0.843222  0.107986         9
#> 6 G00004_drug_004_moa_.. CL00011_cellline_BA_..  0.803893  0.133827         9
#>    fit_source normalization_type
#>   <character>        <character>
#> 1        demo                 GR
#> 2        demo                 RV
#> 3        demo                 GR
#> 4        demo                 RV
#> 5        demo                 GR
#> 6        demo                 RV

One row per (drug × cell line × normalization_type) triplet. The native "Metrics" assay is untouched — because we chose a custom assay name.

Writing to the Metrics assay

You can also write directly to "Metrics" — for example when replacing the standard gDR Hill fit with your own model:

# SE straight from the pipeline — already has a native "Metrics" assay
# (fit_source = "gDR")
"Metrics" %in% assayNames(sa_se)
#> [1] TRUE

# Apply a custom fit to the same assay, coexisting alongside gDR rows
custom_hill <- apply_fit(
  sa_se,
  fit_fn = fit_drug_response_metrics,
  data_type = "single-agent",
  output_assay = "Metrics",
  fit_source = "custom_hill"   # distinct key keeps native "gDR" rows intact
)

metrics_df <- BumpyMatrix::unsplitAsDataFrame(
  assay(custom_hill, "Metrics"),
  row.field = "row", column.field = "column"
)
unique(metrics_df$fit_source)   # both "gDR" and "custom_hill"
#> [1] "gDR"         "custom_hill"

With merge = "merge" (default) the upsert key is fit_source + normalization_type, so native rows (fit_source = "gDR") are preserved. Use merge = "replace" only if you intend to overwrite the whole assay.

Idempotent merge

Calling again with the same fit_source replaces those rows rather than appending:

n_before <- NROW(BumpyMatrix::unsplitAsDataFrame(
  assay(sa_out, "custom_summary"),
  row.field = "row", column.field = "column"
))

# Call again — same fit_source, same data
sa_out2 <- apply_fit(
  sa_out,
  fit_fn = summary_fn,
  data_type = "single-agent",
  output_assay = "custom_summary",
  fit_source = "demo"
)
n_after <- NROW(BumpyMatrix::unsplitAsDataFrame(
  assay(sa_out2, "custom_summary"),
  row.field = "row", column.field = "column"
))

stopifnot(n_before == n_after)   # no duplicate rows
message("Row count before: ", n_before, " — after: ", n_after, " (no change)")

Coexisting fit sources

Different fit_source values live side by side in the same assay:

extra_fn <- function(avg_dt) list(x_max = max(avg_dt$x, na.rm = TRUE))

sa_two <- sa_out |>
  apply_fit(extra_fn, "single-agent",
                   output_assay = "custom_summary",
                   fit_source = "extremes")

sources <- unique(BumpyMatrix::unsplitAsDataFrame(
  assay(sa_two, "custom_summary"),
  row.field = "row", column.field = "column"
)$fit_source)
message("fit_source values in assay: ", paste(sources, collapse = ", "))

Reference Hill fit: fit_drug_response_metrics()

fit_drug_response_metrics() is a reference single-agent fit function that replicates the standard fit_SE() / logisticFit() output exactly.

Property	Value
Model	3-parameter log-logistic (`drc::LL.3u`, `x_0` fixed at 1)
Equivalent to	`fit_SE()` / `gDRutils::logisticFit()`
`x_mean`	Predicted from fitted curve (matches `logisticFit` behaviour)
`fit_type` value	`"DRC3pHillFitModelFixS0"` or `"DRCConstantFitResult"`
Output columns	Full `Metrics` assay schema including `p_value`, `rss`, `x_AOC_range`, `x_max`, `x_sd_avg`

Numerical equivalence with fit_SE()

fit_drug_response_metrics() is numerically identical to fit_SE(). The fit is deterministic — drc::drm does not use random number generation, so set.seed() has no effect. Differences between the two are purely algorithmic:

Same drc::LL.3u model, x_0 = 1, identical priors and concentration bounds
Same x_mean: predicted from the fitted curve over the observed concentration range
Same pcutoff = 0.05 fallback: when the F-test gives p_value ≥ pcutoff, the result is replaced by a flat DRCConstantFitResult (same as logisticFit())
Same n_point_cutoff = 4: fewer unique concentrations → constant fit without attempting the sigmoidal model
All output columns present in the native "Metrics" assay

These parameters can be overridden if your data requires it:

# Loosen the significance threshold or force the sigmoidal fit
fit_drug_response_metrics(avg_dt, pcutoff = 0.1)
fit_drug_response_metrics(avg_dt, force_fit = TRUE)
# Different range for x_AOC_range computation
fit_drug_response_metrics(avg_dt, range_conc = c(1e-3, 10))

hill_out <- apply_fit(
  sa_se,
  fit_fn = fit_drug_response_metrics,
  data_type = "single-agent",
  output_assay = "custom_hill",
  fit_source = "hill_ref"
)

hill_df <- BumpyMatrix::unsplitAsDataFrame(
  assay(hill_out, "custom_hill"),
  row.field = "row", column.field = "column"
)
head(hill_df[, c("row", "column", "normalization_type",
                  "ec50", "xc50", "h", "r2", "fit_type")])
#> DataFrame with 6 rows and 8 columns
#>                      row                 column normalization_type       ec50
#>              <character>            <character>        <character>  <numeric>
#> 1 G00002_drug_002_moa_.. CL00011_cellline_BA_..                 GR 0.00451256
#> 2 G00002_drug_002_moa_.. CL00011_cellline_BA_..                 RV 0.00379589
#> 3 G00003_drug_003_moa_.. CL00011_cellline_BA_..                 GR 0.00471408
#> 4 G00003_drug_003_moa_.. CL00011_cellline_BA_..                 RV 0.00396299
#> 5 G00004_drug_004_moa_.. CL00011_cellline_BA_..                 GR 0.01715280
#> 6 G00004_drug_004_moa_.. CL00011_cellline_BA_..                 RV 0.01654638
#>         xc50         h        r2               fit_type
#>    <numeric> <numeric> <numeric>            <character>
#> 1 0.00895622   1.91190  0.988441 DRC3pHillFitModelFixS0
#> 2 0.00725317   1.82703  0.992470 DRC3pHillFitModelFixS0
#> 3 0.00659147   2.27009  0.990529 DRC3pHillFitModelFixS0
#> 4 0.00565558   2.34960  0.995851 DRC3pHillFitModelFixS0
#> 5        Inf   3.14371  0.991428 DRC3pHillFitModelFixS0
#> 6        Inf   3.15994  0.991887 DRC3pHillFitModelFixS0

ec50 is the raw model parameter (concentration at half-maximal effect); xc50 is the capped version (capped at capping_fold × max(Concentration)) as reported in the native "Metrics" assay.

summary_fn: cell-level aggregation

An optional summary_fn is called once per (drug × cell line) cell on all rows produced by fit_fn for that cell (one row per normalization type). It is the right place for metrics that aggregate across normalization types — for example, whether the result is synergistic across both GR and RV.

# Aggregate: mean xc50 and flag whether both norm types fitted successfully
hill_summary_fn <- function(fit_dt) {
  list(
    mean_xc50 = mean(fit_dt$xc50, na.rm = TRUE),
    mean_r2 = mean(fit_dt$r2, na.rm = TRUE),
    all_converged = all(fit_dt$fit_type == "DRC3pHillFitModelFixS0", na.rm = TRUE)
  )
}

hill_with_summary <- apply_fit(
  sa_se,
  fit_fn = fit_drug_response_metrics,
  data_type = "single-agent",
  output_assay = "custom_hill",
  summary_fn = hill_summary_fn,
  summary_assay = "custom_hill_summary",
  fit_source = "hill_ref"
)
assayNames(hill_with_summary)
#> [1] "RawTreated"          "Controls"            "Normalized"         
#> [4] "Averaged"            "Metrics"             "custom_hill"        
#> [7] "custom_hill_summary"

sumdf <- BumpyMatrix::unsplitAsDataFrame(
  assay(hill_with_summary, "custom_hill_summary"),
  row.field = "row", column.field = "column"
)
head(sumdf[, c("row", "column", "mean_xc50", "mean_r2", "all_converged")])
#> DataFrame with 6 rows and 5 columns
#>                      row                 column  mean_xc50   mean_r2
#>              <character>            <character>  <numeric> <numeric>
#> 1 G00002_drug_002_moa_.. CL00011_cellline_BA_.. 0.00810470  0.990455
#> 2 G00003_drug_003_moa_.. CL00011_cellline_BA_.. 0.00612352  0.993190
#> 3 G00004_drug_004_moa_.. CL00011_cellline_BA_..        Inf  0.991658
#> 4 G00005_drug_005_moa_.. CL00011_cellline_BA_..        Inf  0.977308
#> 5 G00006_drug_006_moa_.. CL00011_cellline_BA_..        Inf  0.997368
#> 6 G00007_drug_007_moa_.. CL00011_cellline_BA_..        Inf  0.850211
#>   all_converged
#>       <logical>
#> 1          TRUE
#> 2          TRUE
#> 3          TRUE
#> 4          TRUE
#> 5          TRUE
#> 6          TRUE

One row per (drug × cell line) — regardless of how many normalization types were fitted.

Combination data: synergy scores

What fit_SE.combinations() does internally

Before describing the new extension API, it helps to understand what the standard fit_SE.combinations() function does. It is a single loop over each (drug-combo × cell-line) pair that executes five sequential steps, each writing to a separate assay:

Step	Key functions	Output assay	Description
1	`fit_combo_cotreatments()`, `fit_combo_codilutions()`	`Metrics`	Fit SA dose-response curves at each co-treatment concentration; produces `ec50`, `h`, `x_inf`, `x_0` per SA series
2	`map_ids_to_fits()`	(internal)	Predict smooth single-agent responses at every combo concentration using the SA fits from step 1; average col/row/codilution predictions
3	`calculate_HSA()`, `calculate_Bliss()`, `calculate_excess()`	`excess`	Compute expected response (HSA = min of SAs; Bliss = product/GR formula); compute per-point excess = expected − observed
4	`calculate_Loewe()`	`isobolograms`, `all_iso_points`	Compute combination index (CI) via isobologram analysis; CI < 1 = synergy
5	`calculate_score()`	`scores`	Reduce per-point excess to a scalar: mean of top-10-percentile values → `bliss_score`, `hsa_score`, `CIScore_50`, `CIScore_80`

The key insight: step 2 (smooth) depends on step 1 (SA fits), and steps 3–5 all depend on step 2. The boundaries between steps are clean, which means they can be extracted as independent public functions.

All five steps are now available as independent public functions:

Step	Public function	Output assay
1+2	`apply_combo_sa_fits()`	`Metrics`
3	`apply_combo_excess()`	`excess`
4	`apply_combo_isobolograms()`	`isobolograms`, `all_iso_points`
5	`apply_combo_scores()`	`scores`

fit_SE.combinations() is now a thin wrapper that calls these four functions in sequence. You can run any subset, swap out a step, or pass a pre-computed assay from a previous step to the next one.

Running all steps individually

# Start from an Averaged-only SE
combo_mae <- gDRutils::get_synthetic_data("finalMAE_combo_matrix_small")
combo_name <- gDRutils::get_supported_experiments("combo")
se <- combo_mae[[combo_name]]
SummarizedExperiment::assays(se) <- SummarizedExperiment::assays(se)["Averaged"]

# Step 1+2: SA fits → Metrics
se <- apply_combo_sa_fits(se)
# Step 3: smooth → excess
se <- apply_combo_excess(se)
# Step 4: Loewe CI → isobolograms
se <- apply_combo_isobolograms(se)
# Step 5: scores
se <- apply_combo_scores(se, excess_assay = "excess")
assayNames(se)
#> [1] "Averaged"       "Metrics"        "excess"         "isobolograms"  
#> [5] "all_iso_points" "scores"

This produces identical results to fit_SE.combinations().

Skipping a step

If you only need SA fits and scores (skipping isobolograms):

se2 <- combo_mae[[combo_name]]
SummarizedExperiment::assays(se2) <- SummarizedExperiment::assays(se2)["Averaged"]

se2 <- apply_combo_sa_fits(se2)
se2 <- apply_combo_excess(se2)
se2 <- apply_combo_scores(se2, excess_assay = "excess")
assayNames(se2)  # no "isobolograms" or "all_iso_points"
#> [1] "Averaged" "Metrics"  "excess"   "scores"

Two approaches for combination scoring in the new API

Two approaches are available, depending on whether you have SA fits or only raw averaged data.

Approach 1: apply_combo_scores() — replicate fit_SE.combinations exactly

apply_combo_scores() is a high-level function that reproduces the Bliss and HSA scoring logic of fit_SE.combinations() exactly, using fitted SA curves from the Metrics assay to generate smooth single-agent predictions.

Property	Value
Requires	`Averaged` + `Metrics` assay (with `dilution_drug`, `ec50`, `h`, …)
Equivalent to	`fit_SE.combinations()` Bliss and HSA scores
Accuracy	cor > 0.998 with `fit_SE.combinations` on real data
Use when	Replacing or extending `fit_SE.combinations` for standard data

# Use the small synthetic combo dataset which has both Averaged and Metrics
combo_mae <- gDRutils::get_synthetic_data("finalMAE_combo_matrix_small")
combo_name <- gDRutils::get_supported_experiments("combo")
combo_se_full <- combo_mae[[combo_name]]

# combo_se_full already has Metrics from fit_SE.combinations
combo_scored <- apply_combo_scores(combo_se_full)
assayNames(combo_scored)
#> [1] "RawTreated"     "Controls"       "Normalized"     "Averaged"      
#> [5] "excess"         "all_iso_points" "isobolograms"   "scores"        
#> [9] "Metrics"

scores_df <- BumpyMatrix::unsplitAsDataFrame(
  assay(combo_scored, "scores"),
  row.field = "row", column.field = "column"
)
scores_df[, c("row", "column", "normalization_type", "bliss_score", "hsa_score")]
#> DataFrame with 24 rows and 5 columns
#>                        row                 column normalization_type
#>                <character>            <character>        <character>
#> 1   G00004_drug_004_moa_.. CL00016_cellline_GB_..                 GR
#> 2   G00004_drug_004_moa_.. CL00016_cellline_GB_..                 RV
#> 3   G00004_drug_004_moa_.. CL00016_cellline_GB_..                 GR
#> 4   G00004_drug_004_moa_.. CL00016_cellline_GB_..                 RV
#> 5   G00005_drug_005_moa_.. CL00016_cellline_GB_..                 GR
#> ...                    ...                    ...                ...
#> 20  G00005_drug_005_moa_.. CL00017_cellline_HB_..                 RV
#> 21  G00006_drug_006_moa_.. CL00017_cellline_HB_..                 GR
#> 22  G00006_drug_006_moa_.. CL00017_cellline_HB_..                 RV
#> 23  G00006_drug_006_moa_.. CL00017_cellline_HB_..                 GR
#> 24  G00006_drug_006_moa_.. CL00017_cellline_HB_..                 RV
#>     bliss_score   hsa_score
#>       <numeric>   <numeric>
#> 1   0.008339057 0.008339057
#> 2   0.002358668 0.002358668
#> 3   0.000674575 0.000674575
#> 4   0.000298478 0.000298478
#> 5   0.009101157 0.009101157
#> ...         ...         ...
#> 20  0.000151856 0.000151856
#> 21  0.031880296 0.031880296
#> 22  0.012769276 0.012769276
#> 23  0.024087057 0.024087057
#> 24  0.005341883 0.005341883

bliss_score > 0 and hsa_score > 0 indicate synergy.

Approach 2: bliss_fit_fn / hss_fit_fn — simplified, no SA fits needed

bliss_fit_fn() and hss_fit_fn() are lower-level fit functions for apply_fit() that compute synergy scores directly from the raw Averaged data, without requiring prior SA curve fits.

Property	Value
Requires	`Averaged` assay only
SA response	Raw single-agent edge points (no curve smoothing)
Use when	Prototyping, custom models, or when SA fits are unavailable

# Build a minimal synthetic combination SE (Averaged only, no Metrics needed)
combo_dt <- data.table::CJ(
  row = c("DrugA", "DrugB"),
  column = "CellLine1",
  normalization_type = c("GR", "RV"),
  Concentration = c(0, 0.1, 1.0),
  Concentration_2 = c(0, 0.1, 1.0)
)
set.seed(42L)
combo_dt[, x := pmax(0.05,
  1 - 0.3 * Concentration / (Concentration + 0.5) -
      0.2 * Concentration_2 / (Concentration_2 + 0.5) +
      rnorm(.N, 0, 0.03))]

data_cols <- setdiff(names(combo_dt), c("row", "column"))
combo_bumpy <- BumpyMatrix::splitAsBumpyMatrix(
  combo_dt[, data_cols, with = FALSE],
  row = combo_dt$row, col = combo_dt$column
)
combo_se <- SummarizedExperiment(assays = list(Averaged = combo_bumpy))

bliss_out <- apply_fit(
  combo_se,
  fit_fn = bliss_fit_fn,
  data_type = "combination",
  output_assay = "custom_bliss",
  fit_source = "bliss"
)

bliss_df <- BumpyMatrix::unsplitAsDataFrame(
  assay(bliss_out, "custom_bliss"),
  row.field = "row", column.field = "column"
)
bliss_df[, c("row", "column", "normalization_type",
             "bliss_score", "bliss_excess_mean", "n_combo_points")]
#> DataFrame with 4 rows and 6 columns
#>           row      column normalization_type bliss_score bliss_excess_mean
#>   <character> <character>        <character>   <numeric>         <numeric>
#> 1       DrugA   CellLine1                 GR   0.0404576        0.02194040
#> 2       DrugA   CellLine1                 RV   0.1790833        0.07174050
#> 3       DrugB   CellLine1                 GR   0.1083000        0.02410753
#> 4       DrugB   CellLine1                 RV   0.0473357        0.00199227
#>   n_combo_points
#>        <integer>
#> 1              4
#> 2              4
#> 3              4
#> 4              4

hss_out <- apply_fit(
  combo_se,
  fit_fn = hss_fit_fn,
  data_type = "combination",
  output_assay = "custom_hss",
  fit_source = "hss"
)

hss_df <- BumpyMatrix::unsplitAsDataFrame(
  assay(hss_out, "custom_hss"),
  row.field = "row", column.field = "column"
)
hss_df[, c("row", "column", "normalization_type",
           "hss_score", "hss_excess_mean")]
#> DataFrame with 4 rows and 5 columns
#>           row      column normalization_type hss_score hss_excess_mean
#>   <character> <character>        <character> <numeric>       <numeric>
#> 1       DrugA   CellLine1                 GR  0.118126       0.0711634
#> 2       DrugA   CellLine1                 RV  0.232106       0.0971858
#> 3       DrugB   CellLine1                 GR  0.197907       0.0726057
#> 4       DrugB   CellLine1                 RV  0.166576       0.0480229

Efficient multi-fit: apply_fits()

When several fit functions operate on the same input assay, use apply_fits() to traverse each BumpyMatrix cell once and apply all functions in that single pass.

combo_multi <- apply_fits(
  combo_se,
  fit_fns = list(
    custom_bliss = bliss_fit_fn,
    custom_hss = hss_fit_fn
  ),
  data_type = "combination",
  fit_source = "synergy_panel"
)
assayNames(combo_multi)
#> [1] "Averaged"     "custom_bliss" "custom_hss"

Both assays are written in one traversal — equivalent to chaining two apply_fit() calls but without the overhead of a second unsplit + iteration.

Shared pre-computation pattern

When two metrics share an expensive intermediate (e.g. fitted single-agent curves), a single fit function can return a named list of named lists to populate multiple assays from one computation:

# Each top-level name maps to an output assay; inner lists are the rows
bliss_and_hss_combined <- function(dt) {
  # Expensive step done ONCE per cell
  sa1 <- dt[dt$Concentration_2 == 0 & dt$Concentration > 0, ]
  sa2 <- dt[dt$Concentration == 0   & dt$Concentration_2 > 0, ]

  list(
    custom_bliss = list(bliss_score = mean(sa1$x) - mean(sa2$x)), # simplified
    custom_hss = list(hss_score = min(c(sa1$x, sa2$x)))
  )
}

apply_fits(
  combo_se,
  fit_fns = list(custom_bliss = bliss_and_hss_combined,
                    custom_hss = bliss_and_hss_combined),
  data_type = "combination",
  fit_source = "shared"
)

Chaining with pipe

The functions are pipe-friendly — each call returns the updated SE:

result_se <- combo_se |>
  apply_fit(bliss_fit_fn, "combination",
                   output_assay = "custom_bliss",
                   fit_source = "bliss") |>
  apply_fit(hss_fit_fn, "combination",
                   output_assay = "custom_hss",
                   fit_source = "hss") |>
  apply_fit(
    function(dt) list(n_obs = NROW(dt)),
    "combination",
    output_assay = "combo_diagnostics",
    fit_source = "qc"
  )

assayNames(result_se)
#> [1] "Averaged"          "custom_bliss"      "custom_hss"       
#> [4] "combo_diagnostics"

Error handling

By default, a failed cell emits a warning and is skipped (on_error = "warn"). Use on_error = "stop" to halt immediately and propagate the error — useful when debugging a new fit function.

buggy_fn <- function(dt) {
  if (dt$normalization_type[1] == "GR") stop("GR not supported")
  list(x_rv = mean(dt$x, na.rm = TRUE))
}

se_partial <- withCallingHandlers(
  apply_fit(
    sa_se, buggy_fn, "single-agent",
    output_assay = "rv_only",
    fit_source = "rv_fn",
    on_error = "warn"
  ),
  warning = function(w) {
    message("[caught] ", conditionMessage(w))
    invokeRestart("muffleWarning")
  }
)
# Only RV rows are written; GR cells were skipped with a warning
rv_df <- BumpyMatrix::unsplitAsDataFrame(
  assay(se_partial, "rv_only"),
  row.field = "row", column.field = "column"
)
unique(rv_df$normalization_type)
#> [1] "RV"

Time-course data: `fit_SE.timecourse()`

Live-cell imaging experiments (e.g. Incucyte) produce per-well cell counts at multiple time points rather than endpoint readouts. The standard fit_SE() pipeline is not appropriate for this layout. fit_SE.timecourse() provides a dedicated two-stage pipeline:

Stage	What it does	Output
1	`lm(LogFoldChange ~ Duration)` per well per time window → mean growth rate (doublings/day) across replicates	`GrowthRates` assay
2	Logistic Hill fit on `NormalizedGrowthRate ~ log10(Concentration)` per `(Drug, CellLine, period)`	`Metrics` assay

Inputs

fit_SE.timecourse() expects a SummarizedExperiment with a LogFoldChange assay — the output of normalize_SE(se, data_type = "time-course"):

# periods: named list of half-open c(start_hour, end_hour) windows,
# i.e. a timepoint is used when start_hour <= Duration < end_hour
periods <- list(
  early = c(0, 48),
  late = c(48, 96)
)

# normalization_map: per-period DMSO reference; must name every period
#   "None"  → keep raw growth rate (doublings/day)
#   "early" → divide by DMSO rate from the "early" period
norm_map <- c(early = "None", late = "early")

se_fit <- fit_SE.timecourse(
  se_tc,                          # SE with LogFoldChange assay
  periods = periods,
  normalization_map = norm_map
)

Output structure

The result is a SummarizedExperiment where columns = period names:

# One column per period
colnames(se_fit)   # c("early", "late")

# Metrics assay: fitted parameters per (Drug × CellLine × period)
metrics_dt <- gDRutils::convert_se_assay_to_dt(se_fit, "Metrics")
head(metrics_dt[, c("DrugName", "CellLineName", "ec50", "xc50", "r2",
                     "x_inf", "x_0", "N_conc")])

Each row is one dose-response curve for a (Drug, CellLine, period) combination, with the same fit columns as fit_SE(): ec50, xc50, h, r2, x_mean, x_AOC, N_conc, fit_type, etc.

Custom fitting for time-course data

Both stages are pluggable, so the custom fit interface covers time-course data the same way it covers single-agent and combination screens.

Stage 2 takes a fit_fn with the usual contract — it receives one normalization_type = "NGR" slice of the GrowthRates assay, with NormalizedGrowthRate as the response and Concentration as the x axis:

auc_fn <- function(dt) {
  list(auc_ngr = mean(dt$NormalizedGrowthRate, na.rm = TRUE))
}

se_fit <- fit_SE.timecourse(
  se_tc,
  periods = periods,
  normalization_map = norm_map,
  fit_fn = auc_fn                 # fit_source defaults to "custom"
)

Stage 1 takes a rate_fn, which receives the timepoints of a single well within a single window and returns one growth rate in doublings per day. Replace it when the default lm(LogFoldChange ~ Duration) slope is not the right estimator — for example for screens with only two usable timepoints, or when the growth phase has to be selected per well:

endpoint_rate_fn <- function(dt) {
  ord <- order(dt$Duration)
  24 * diff(range(dt$LogFoldChange[ord])) / diff(range(dt$Duration[ord]))
}

se_fit <- fit_SE.timecourse(
  se_tc,
  periods = periods,
  normalization_map = norm_map,
  rate_fn = endpoint_rate_fn
)

The stages can also be run separately, which is the same two-step workflow as average_SE() → apply_fit() for single-agent data. compute_growth_rates() returns the growth rate table — one row per (CellLine, Drug, Concentration, partner slots, period), controls included, with GrowthRate, sd_GrowthRate and rate_0 — and growth_rates_to_se() turns it into the stage 2 input:

growth_dt <- compute_growth_rates(se_tc, periods, norm_map, rate_fn = endpoint_rate_fn)

se_fit <- apply_fit(
  growth_rates_to_se(growth_dt),
  fit_fn = auc_fn,
  data_type = "time-course-metrics", # profile supplies input_assay = "GrowthRates"
  output_assay = "Metrics",
  merge = "replace",
  fit_source = "custom"
)

Note that fit_SE.timecourse() keeps the GrowthRates assay in its output, so an already-fitted object can be re-fitted with a different fit_fn without recomputing stage 1.

Key design notes

normalization_type = "NGR" (NormalizedGrowthRate) — distinct from "GR" (GR value) and "RV" (Relative Viability) used in single-agent screens. Time-course data has no RV/GR split.
Half-open windows [start, end) — contiguous periods partition the time-course: the shared boundary belongs to the later window and t = 0 is kept.
All partner-drug slots in the row key — each combination arm at each partner concentration gets its own dose-response curve.
untreated_tag from gDRutils — DMSO detection uses get_env_identifiers("untreated_tag") rather than hardcoded strings.
Two profiles, one apply_fit() pass — "time-course" supplies the default input assay of stage 1, "time-course-metrics" drives the stage 2 slicing. Stage 1 aggregates replicates and normalises against control wells that live in other BumpyMatrix cells, so it cannot be expressed as an apply_fit() pass; that is why it is pluggable through rate_fn instead. See GDR-3506 for the planned multi-stage profile model.

Quick reference

Function signatures

# Single fit → one output assay
apply_fit(
  se,
  fit_fn,
  data_type = "single-agent", # or "combination", "time-course-metrics"
  slicing_cols = NULL, # NULL → data_type default
  slicing_values = NULL, # NULL → all unique values found
  input_assay = NULL, # NULL → data_type default ("Averaged")
  output_assay, # REQUIRED — your assay name
  summary_fn = NULL, # optional cell-level aggregator
  summary_assay = NULL,
  merge = "merge", # or "replace"
  on_error = "warn", # or "stop"
  fit_source                         # REQUIRED — upsert key tag
)

# Multiple fits → one BumpyMatrix pass
apply_fits(
  se,
  fit_fns, # named list: name = output assay, value = fit function
  data_type = "single-agent",
  fit_source,
  ...
)

fit_fn contract

Input	`data.table` — one BumpyMatrix cell, filtered to one `slicing_cols` value
Output	Named list → one row in `output_assay`; or named list of named lists for multi-assay pattern
fit_source	Stamped by the generic layer — do not include in the return value
slicing columns	The caller’s value (e.g. `normalization_type`) IS available in the data.table — no need to filter again

summary_fn contract

Input	`data.table` — all rows written by `fit_fn` for one (row × column) cell
Output	Named list → one row in `summary_assay`
When to use	Aggregated metrics that span normalization types (e.g. “converged in at least one type?”)

Reference implementations

Single-agent fit functions (for `apply_fit()`)

Function	Model	Equivalent to	Key output columns
`fit_drug_response_metrics()`	3p LL.3u, `x_0 = 1`	`fit_SE()` / `logisticFit()`	`ec50`, `xc50`, `h`, `r2`, `x_mean`, `x_AOC`, `fit_type = "DRC3pHillFitModelFixS0"`
`fit_drug_response_metrics_4p()`	4p LL.4, `x_0` free	— (extended variant)	`ec50`, `xc50`, `h`, `r2`, `x_0`, `x_mean`, `fit_type = "DRC4pHillFitModel"`

Combination scoring

Composable step functions

Function	Step	Output assay	Description
`apply_combo_sa_fits(se)`	1+2	`Metrics`	SA curve fits + smooth predictions
`apply_combo_excess(se)`	3	`excess`	HSA/Bliss excess per combo point
`apply_combo_isobolograms(se)`	4	`isobolograms`, `all_iso_points`	Loewe CI isobologram
`apply_combo_scores(se)`	5	`scores`	Bliss/HSA scalar scores

fit_SE.combinations() is a thin wrapper that calls all four in order.

Lower-level scoring functions (for `apply_fit()`)

Function	Level	Requires	Equivalent to	Key output columns
`apply_combo_scores(se)`	SE-level (recommended)	`Averaged` + `Metrics`	`fit_SE.combinations()` Bliss & HSA	`bliss_score`, `hsa_score`
`bliss_fit_fn()`	triplet `fit_fn`	`Averaged` only	— (simplified, no SA fits)	`bliss_score`, `bliss_excess_mean`, `n_combo_points`
`hss_fit_fn()`	triplet `fit_fn`	`Averaged` only	— (simplified, no SA fits)	`hss_score`, `hss_excess_mean`, `n_combo_points`

When to use which: - Use apply_combo_scores() when you have a fully fitted SE (after fit_SE.combinations() or apply_fit_to_se()) and want scores numerically consistent with the standard gDR pipeline. - Use bliss_fit_fn() / hss_fit_fn() when prototyping a new scoring approach, when SA fits are unavailable, or when embedding score computation inside a larger custom fit_fn.

SessionInfo

sessionInfo()
#> R version 4.6.1 (2026-06-24)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 26.04.1 LTS
#> 
#> Matrix products: default
#> BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
#> LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.32.so;  LAPACK version 3.12.0
#> 
#> locale:
#>  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
#>  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
#>  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
#>  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
#>  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
#> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
#> 
#> time zone: Etc/UTC
#> tzcode source: system (glibc)
#> 
#> attached base packages:
#> [1] stats4    stats     graphics  grDevices utils     datasets  methods  
#> [8] base     
#> 
#> other attached packages:
#>  [1] data.table_1.18.6.1         BumpyMatrix_1.21.0         
#>  [3] SummarizedExperiment_1.43.0 Biobase_2.73.2             
#>  [5] GenomicRanges_1.65.4        Seqinfo_1.3.2              
#>  [7] IRanges_2.47.5              S4Vectors_0.51.10          
#>  [9] BiocGenerics_0.59.12        generics_0.1.4             
#> [11] MatrixGenerics_1.25.0       matrixStats_1.5.0          
#> [13] gDRutils_1.11.17            gDRtestData_1.11.10        
#> [15] gDRcore_1.11.18             BiocStyle_2.41.0           
#> 
#> loaded via a namespace (and not attached):
#>  [1] xfun_0.61                   bslib_0.12.0               
#>  [3] lattice_0.23-1              tools_4.6.1                
#>  [5] parallel_4.6.1              sandwich_3.1-3             
#>  [7] drc_4.0-0                   Matrix_1.7-6               
#>  [9] checkmate_2.3.4             RColorBrewer_1.1-3         
#> [11] RcppParallel_6.2.1          lifecycle_1.0.5            
#> [13] farver_2.1.2                compiler_4.6.1             
#> [15] codetools_0.2-20            carData_3.0-6              
#> [17] htmltools_0.5.9             sys_3.4.3                  
#> [19] buildtools_1.0.0            sass_0.4.10                
#> [21] yaml_2.3.12                 Formula_1.2-6              
#> [23] car_3.1-5                   jquerylib_0.1.4            
#> [25] MASS_7.3-66                 BiocParallel_1.47.0        
#> [27] DelayedArray_0.39.8         cachem_1.1.0               
#> [29] abind_1.4-8                 multcomp_1.4-32            
#> [31] qs2_0.3.1                   gtools_3.9.5               
#> [33] digest_0.6.39               mvtnorm_1.4-2              
#> [35] maketools_1.3.2             splines_4.6.1              
#> [37] fastmap_1.2.0               grid_4.6.1                 
#> [39] cli_3.6.6                   SparseArray_1.13.4         
#> [41] S4Arrays_1.13.2             survival_3.8-12            
#> [43] TH.data_1.1-5               scales_1.4.0               
#> [45] backports_1.5.1             plotrix_3.8-14             
#> [47] lambda.r_1.2.4              rmarkdown_2.32             
#> [49] XVector_0.53.0              otel_0.2.0                 
#> [51] futile.logger_1.4.9         zoo_1.9-1                  
#> [53] stringfish_0.19.2           evaluate_1.0.5             
#> [55] knitr_1.52                  MultiAssayExperiment_1.39.1
#> [57] futile.options_1.0.1        rlang_1.3.0                
#> [59] Rcpp_1.1.2                  glue_1.8.1                 
#> [61] formatR_1.14                BiocManager_1.30.27        
#> [63] jsonlite_2.0.0              R6_2.6.1

Custom fitting with apply_fit()