This report summarizes the data available in the curated CTRPv2 PharmacoSet, with emphasis on sample coverage, molecular-profile distributions, treatment response, and agreement between published CTRP metrics and values recomputed with PharmacoGx.

Load the PharmacoSet

pset <- readRDS(params$pset_path)
pset
## <PharmacoSet>
## Name: CTRPv2 
## Date Created: Thu Jul 30 13:05:24 2026 
## Number of samples:  1107 
## Molecular profiles: <MultiAssayExperiment> 
##    ExperimentList class object of length 2: 
##     [1] microarray.gene_expression_log2 : SummarizedExperiment with 17970 rows and 823 columns 
##     [2] cnv.gene_log2 : SummarizedExperiment with 22598 rows and 823 columns 
## Treatment response: <TreatmentResponseExperiment> 
##    dim:  395263 1107 
##    assays(2): profiles sensitivity 
##    rownames(395263): CTRPv2.0_EXP_100_CPD_100490 CTRPv2.0_EXP_100_CPD_102193 ... CTRPv2.0_EXP_9_CPD_705300 CTRPv2.0_EXP_9_CPD_96171 
##    rowData(12): response_id experiment_id master_cpd_id ... cells_per_well growth_mode snp_fp_status 
##    colnames(1107): 1 3 4 ... 155826 155827 155841 
##    colData(1): master_ccl_id 
##    metadata(7): data_source ... recomputed_profiles_in_profiles_assay created sessionInfo raw_dose_response_rows ... recomputed_profiles_in_profiles_assay created sessionInfo aggregated_dose_response_rows ... recomputed_profiles_in_profiles_assay created sessionInfo

Sample and Assay Coverage

The PharmacoSet combines sample metadata with treatment-response, microarray, and copy-number data. Figure 1 shows both the size of each data layer and the exact intersections among them.

samples <- as.data.table(PharmacoGx::sampleInfo(pset))
treatments <- as.data.table(PharmacoGx::treatmentInfo(pset))
molecular_profiles <- PharmacoGx::molecularProfiles(pset)
treatment_response <- PharmacoGx::treatmentResponse(pset)

microarray_profile <- MultiAssayExperiment::experiments(molecular_profiles)[[
  "microarray.gene_expression_log2"
]]
cnv_profile <- MultiAssayExperiment::experiments(molecular_profiles)[[
  "cnv.gene_log2"
]]
profiles <- as.data.table(CoreGx::assay(treatment_response, "profiles"))

sample_availability <- data.table(
  master_ccl_id = as.character(samples$master_ccl_id)
)
sample_availability[,
  treatment_response := master_ccl_id %in%
    as.character(profiles$master_ccl_id)
]
sample_availability[,
  microarray := master_ccl_id %in%
    as.character(colnames(microarray_profile))
]
sample_availability[,
  cnv := master_ccl_id %in%
    as.character(colnames(cnv_profile))
]
sample_availability[, metadata := TRUE]

upset_input <- sample_availability[,
  .(
    Metadata = as.integer(metadata),
    `Treatment response` = as.integer(treatment_response),
    Microarray = as.integer(microarray),
    CNV = as.integer(cnv)
  )
]

UpSetR::upset(
  as.data.frame(upset_input),
  sets = c(
    "Metadata",
    "Treatment response",
    "Microarray",
    "CNV"
  ),
  nsets = 4,
  nintersects = NA,
  keep.order = TRUE,
  order.by = "freq",
  decreasing = TRUE,
  show.numbers = "yes",
  number.angles = 0,
  point.size = 3.2,
  line.size = 1,
  main.bar.color = "#3B82A0",
  sets.bar.color = "#6B7280",
  mainbar.y.label = "Samples in intersection",
  sets.x.label = "Samples in data layer",
  text.scale = c(1.25, 1.1, 1.1, 1.05, 1.1, 1.05)
)
Sample intersections across metadata, treatment response, microarray, and CNV layers.

Figure 1: Sample intersections across metadata, treatment response, microarray, and CNV layers.

Of the 1,107 samples in the PharmacoSet, all have metadata, 887 have treatment-response profiles, 823 have microarray data, and 823 have CNV data. 220 samples have metadata without a response or molecular assay. Metadata is a shared layer in the plot, not a separate metadata-only group.

The equal 12-sample intersections are not a metadata alignment artifact. The microarray and CNV assays each contain 823 unique sample identifiers, but only 811 are shared: 12 response-tested samples occur only in microarray and a different 12 occur only in CNV. All 24 identifiers are unique and map directly to the matching PSet sample identifier. The following matrix names those source-specific samples and shows their assay membership.

discordant_samples <- merge(
  sample_availability[
    treatment_response & xor(microarray, cnv)
  ],
  samples[, .(
    master_ccl_id = as.character(master_ccl_id),
    CELL_NAME
  )],
  by = "master_ccl_id",
  all.x = TRUE
)
discordant_samples[,
  coverage_group := fifelse(
    microarray,
    "Only Microarray",
    "Only CNV"
  )
]
discordant_samples[,
  sample_label := paste0(
    CELL_NAME,
    " (",
    master_ccl_id,
    ")"
  )
]
discordant_coverage <- melt(
  discordant_samples,
  id.vars = c("master_ccl_id", "sample_label", "coverage_group"),
  measure.vars = c("treatment_response", "microarray", "cnv"),
  variable.name = "data_layer",
  value.name = "present"
)
discordant_coverage[,
  data_layer := factor(
    data_layer,
    levels = c("treatment_response", "microarray", "cnv"),
    labels = c("Treatment response", "Microarray", "CNV")
  )
]

ggplot(
  discordant_coverage,
  aes(x = data_layer, y = reorder(sample_label, master_ccl_id), fill = present)
) +
  geom_tile(colour = "white", linewidth = 0.8) +
  facet_grid(
    rows = vars(coverage_group),
    scales = "free_y",
    space = "free_y"
  ) +
  scale_fill_manual(
    values = c("FALSE" = "#E5E7EB", "TRUE" = "#3B82A0"),
    labels = c("Absent", "Present"),
    name = NULL
  ) +
  labs(
    title = "Samples with discordant molecular-profile coverage",
    x = NULL,
    y = NULL
  ) +
  theme(
    panel.grid = element_blank(),
    legend.position = "top",
    axis.text.x = element_text(angle = 20, hjust = 1)
  )
Treatment-response samples represented in only one molecular profile.

Figure 2: Treatment-response samples represented in only one molecular profile.

Sample and Treatment Annotation Coverage

AnnotationDB is the primary annotation source. AnnotationGx supplies a fallback when the primary source does not provide the final Cellosaurus accession or PubChem CID needed for cross-dataset use. Figure 3 reports the coverage present in the loaded PSet.

sample_annotation_sources <- samples[,
  .(
    entity = "Samples: Cellosaurus",
    source = fcase(
      !is.na(AnnotationDB.Cellosaurus.Accession) & as.character(AnnotationDB.Cellosaurus.Accession) != "" ,
      "AnnotationDB"                                                                                      ,
      !is.na(AnnotationGx.Cellosaurus.Accession) & as.character(AnnotationGx.Cellosaurus.Accession) != "" ,
      "AnnotationGx fallback"                                                                             ,
      default = "Unresolved"
    )
  )
]

treatment_annotation_sources <- treatments[,
  .(
    entity = "Treatments: PubChem CID",
    source = fcase(
      !is.na(AnnotationDB.PubChem.CID) ,
      "AnnotationDB"                   ,
      !is.na(AnnotationGx.PubChem.CID) ,
      "AnnotationGx fallback"          ,
      default = "Unresolved"
    )
  )
]

annotation_counts <- rbind(
  sample_annotation_sources,
  treatment_annotation_sources
)[,
  .(records = .N),
  by = .(entity, source)
]
annotation_coverage <- merge(
  CJ(
    entity = c("Samples: Cellosaurus", "Treatments: PubChem CID"),
    source = c("AnnotationDB", "AnnotationGx fallback", "Unresolved"),
    unique = TRUE
  ),
  annotation_counts,
  by = c("entity", "source"),
  all.x = TRUE
)
annotation_coverage[is.na(records), records := 0L]
annotation_coverage[, total := sum(records), by = entity]
annotation_coverage[, fraction := records / total]
annotation_coverage[,
  source := factor(
    source,
    levels = c("AnnotationDB", "AnnotationGx fallback", "Unresolved")
  )
]

sample_annotation_caption <- annotation_coverage[
  entity == "Samples: Cellosaurus"
]
treatment_annotation_caption <- annotation_coverage[
  entity == "Treatments: PubChem CID"
]

unresolved_samples <- samples[
  (
    is.na(AnnotationDB.Cellosaurus.Accession) |
      as.character(AnnotationDB.Cellosaurus.Accession) == ""
  ) &
    (
      is.na(AnnotationGx.Cellosaurus.Accession) |
        as.character(AnnotationGx.Cellosaurus.Accession) == ""
    )
]
unresolved_treatments <- treatments[
  is.na(AnnotationDB.PubChem.CID) &
    is.na(AnnotationGx.PubChem.CID)
]
unresolved_treatments[,
  unresolved_reason := fcase(
    !is.na(Compound.Name) &
      grepl(":", Compound.Name, fixed = TRUE),
    "Fixed-ratio mixture",
    !is.na(Compound.Name) &
      grepl("^BRD-[A-Z][0-9]{8}$", Compound.Name),
    "Broad identifier used as name",
    !is.na(source_catalog_id) &
      as.character(source_catalog_id) != "",
    "Source catalog entry without match",
    default = "Named compound without match"
  )
]

ggplot(
  annotation_coverage,
  aes(x = entity, y = records, fill = source)
) +
  geom_col(position = "fill", width = 0.65) +
  geom_text(
    aes(
      label = ifelse(
        fraction >= 0.035,
        paste0(comma(records), " / ", comma(total)),
        ""
      )
    ),
    position = position_fill(vjust = 0.5),
    colour = "white",
    fontface = "bold",
    size = 3.6
  ) +
  scale_y_continuous(labels = percent) +
  scale_fill_manual(
    values = c(
      "AnnotationDB" = "#26734D",
      "AnnotationGx fallback" = "#D99032",
      "Unresolved" = "#9B3A42"
    ),
    drop = FALSE
  ) +
  labs(
    title = "Coverage of external identifiers",
    subtitle = "Unresolved = no identifier from AnnotationDB or AnnotationGx",
    x = NULL,
    y = "Percentage",
    fill = NULL,
    caption = paste0(
      "Samples - AnnotationDB: ",
      comma(sample_annotation_caption[source == "AnnotationDB", records]),
      " / ",
      comma(nrow(samples)),
      "; AnnotationGx fallback: ",
      comma(sample_annotation_caption[
        source == "AnnotationGx fallback",
        records
      ]),
      " / ",
      comma(nrow(samples)),
      "; unresolved: ",
      comma(sample_annotation_caption[source == "Unresolved", records]),
      " / ",
      comma(nrow(samples)),
      ".\nTreatments - AnnotationDB: ",
      comma(treatment_annotation_caption[source == "AnnotationDB", records]),
      " / ",
      comma(nrow(treatments)),
      "; AnnotationGx fallback: ",
      comma(treatment_annotation_caption[
        source == "AnnotationGx fallback",
        records
      ]),
      " / ",
      comma(nrow(treatments)),
      "; unresolved: ",
      comma(treatment_annotation_caption[source == "Unresolved", records]),
      " / ",
      comma(nrow(treatments)),
      ".\nUnresolved samples without a Cellosaurus match: ",
      comma(nrow(unresolved_samples)),
      ".\nUnresolved treatments - fixed-ratio mixtures: ",
      comma(sum(
        unresolved_treatments$unresolved_reason == "Fixed-ratio mixture"
      )),
      ".\nBroad identifiers used as names: ",
      comma(sum(
        unresolved_treatments$unresolved_reason ==
          "Broad identifier used as name"
      )),
      "; source catalog entries: ",
      comma(sum(
        unresolved_treatments$unresolved_reason ==
          "Source catalog entry without match"
      )),
      "; named compounds: ",
      comma(sum(
        unresolved_treatments$unresolved_reason ==
          "Named compound without match"
      )),
      ".\nMixtures generally do not have one appropriate PubChem CID."
    )
  )
Cellosaurus and PubChem CID coverage by annotation source.

Figure 3: Cellosaurus and PubChem CID coverage by annotation source.

Sample Bias Estimation

Sex and parsed population are summarized for the metadata collection and for each data layer. AnnotationDB.Population.Parsed is curated from AnnotationDB.Population; source strings mix population, ancestry, nationality, and combined descriptions. It is not a standardized or self-reported race variable.

bias_samples <- copy(samples)
bias_samples[, master_ccl_id := as.character(master_ccl_id)]
if (!"AnnotationDB.Population.Parsed" %in% names(bias_samples)) {
  stop(
    "Sample metadata is missing AnnotationDB.Population.Parsed.",
    call. = FALSE
  )
}
bias_samples[,
  sex_group := fcase(
    is.na(AnnotationDB.Sex) |
      trimws(as.character(AnnotationDB.Sex)) == "",
    "Missing",
    tolower(trimws(as.character(AnnotationDB.Sex))) == "male",
    "Male",
    tolower(trimws(as.character(AnnotationDB.Sex))) == "female",
    "Female",
    tolower(trimws(as.character(AnnotationDB.Sex))) == "sex unspecified",
    "Sex unspecified",
    default = "Other/unclassified"
  )
]
bias_samples[,
  derived_race := fcoalesce(
    as.character(AnnotationDB.Population.Parsed),
    "Missing"
  )
]

sex_levels <- c(
  "Male",
  "Female",
  "Sex unspecified",
  "Other/unclassified",
  "Missing"
)
race_levels <- c(
  "White",
  "Black/African",
  "East Asian",
  "South Asian",
  "Asian unspecified",
  "Multiracial",
  "Other/unclassified",
  "Missing"
)
bias_samples[,
  sex_group := factor(sex_group, levels = sex_levels)
]
bias_samples[,
  derived_race := factor(derived_race, levels = race_levels)
]

demographic_composition <- rbindlist(list(
  bias_samples[, .(
    dimension = "Sex",
    category = as.character(sex_group)
  )],
  bias_samples[, .(
    dimension = "Derived race",
    category = as.character(derived_race)
  )]
))[
  ,
  .(samples = .N),
  by = .(dimension, category)
]
demographic_composition[,
  dimension := factor(
    dimension,
    levels = c("Sex", "Derived race")
  )
]
demographic_composition[,
  category := factor(
    category,
    levels = unique(c(sex_levels, race_levels))
  )
]

ggplot(
  demographic_composition,
  aes(x = samples, y = category, fill = dimension)
) +
  geom_col(width = 0.72) +
  geom_text(
    aes(label = comma(samples)),
    hjust = -0.12,
    colour = "grey20",
    size = 3.6
  ) +
  facet_grid(
    rows = vars(dimension),
    scales = "free_y",
    space = "free_y"
  ) +
  scale_x_continuous(
    labels = comma,
    expand = expansion(mult = c(0, 0.13))
  ) +
  scale_fill_manual(
    values = c("Sex" = "#4C78A8", "Derived race" = "#9B6A9A"),
    guide = "none"
  ) +
  labs(
    title = "Sample composition",
    x = "Samples",
    y = NULL,
    caption = paste0(
      "Parsed population is curated from AnnotationDB population strings. ",
      "Original strings remain in the PSet. Other/unclassified: ",
      comma(bias_samples[derived_race == "Other/unclassified", .N]),
      "."
    )
  ) +
  theme(
    strip.text.y = element_text(angle = 0),
    panel.grid.major.y = element_blank()
  )
Sex and derived race composition of the PSet sample metadata.

Figure 4: Sex and derived race composition of the PSet sample metadata.

Current metadata contains 601 male and 451 female samples. Largest derived race groups are White (671) and East Asian (302). All nonmissing parsed population values use a documented group or remain visible as Other/unclassified.

Assay coverage within each subgroup is compared with overall layer coverage. Open points mark groups with fewer than 10 samples.

bias_samples <- merge(
  bias_samples,
  sample_availability[, .(
    master_ccl_id,
    treatment_response,
    microarray,
    cnv
  )],
  by = "master_ccl_id",
  all.x = TRUE
)

sex_assay_membership <- melt(
  bias_samples[, .(
    category = as.character(sex_group),
    treatment_response,
    microarray,
    cnv
  )],
  id.vars = "category",
  measure.vars = c("treatment_response", "microarray", "cnv"),
  variable.name = "data_layer",
  value.name = "present"
)
sex_assay_membership[, dimension := "Sex"]

race_assay_membership <- melt(
  bias_samples[, .(
    category = as.character(derived_race),
    treatment_response,
    microarray,
    cnv
  )],
  id.vars = "category",
  measure.vars = c("treatment_response", "microarray", "cnv"),
  variable.name = "data_layer",
  value.name = "present"
)
race_assay_membership[, dimension := "Derived race"]

subgroup_assay_coverage <- rbindlist(list(
  sex_assay_membership,
  race_assay_membership
))[
  ,
  .(
    group_n = .N,
    represented = sum(present),
    coverage = mean(present)
  ),
  by = .(dimension, category, data_layer)
]
overall_assay_coverage <- data.table(
  data_layer = c("treatment_response", "microarray", "cnv"),
  overall_coverage = c(
    mean(sample_availability$treatment_response),
    mean(sample_availability$microarray),
    mean(sample_availability$cnv)
  )
)
subgroup_assay_coverage <- merge(
  subgroup_assay_coverage,
  overall_assay_coverage,
  by = "data_layer",
  all.x = TRUE
)
subgroup_assay_coverage[,
  data_layer := factor(
    data_layer,
    levels = c("treatment_response", "microarray", "cnv"),
    labels = c("Treatment response", "Microarray", "CNV")
  )
]
subgroup_assay_coverage[,
  dimension := factor(
    dimension,
    levels = c("Sex", "Derived race")
  )
]
subgroup_assay_coverage[,
  category_order := fifelse(
    dimension == "Sex",
    match(category, sex_levels),
    match(category, race_levels)
  )
]
subgroup_assay_coverage[,
  group_label := paste0(category, " (n = ", comma(group_n), ")")
]
coverage_label_order <- subgroup_assay_coverage[
  order(dimension, category_order),
  unique(group_label)
]
subgroup_assay_coverage[,
  group_label := factor(
    group_label,
    levels = rev(coverage_label_order)
  )
]
subgroup_assay_coverage[, small_group := group_n < 10]
subgroup_assay_coverage[,
  coverage_difference := coverage - overall_coverage
]

stable_selection_shifts <- subgroup_assay_coverage[group_n >= 10]
largest_stable_shift <- stable_selection_shifts[
  which.max(abs(coverage_difference))
]
selection_shift_statement <- paste0(
  "Largest subgroup shift among groups with at least 10 samples: ",
  largest_stable_shift$category,
  " in ",
  as.character(largest_stable_shift$data_layer),
  " (",
  percent(
    largest_stable_shift$coverage_difference,
    accuracy = 0.1
  ),
  " versus overall coverage)."
)

ggplot(
  subgroup_assay_coverage,
  aes(x = coverage, y = group_label)
) +
  geom_vline(
    aes(xintercept = overall_coverage),
    colour = "grey55",
    linetype = "dashed",
    linewidth = 0.7
  ) +
  geom_segment(
    aes(
      x = overall_coverage,
      xend = coverage,
      yend = group_label
    ),
    colour = "grey75",
    linewidth = 0.7
  ) +
  geom_point(
    aes(shape = small_group),
    colour = "#2D6584",
    fill = "#2D6584",
    size = 3
  ) +
  facet_grid(
    rows = vars(dimension),
    cols = vars(data_layer),
    scales = "free_y",
    space = "free_y"
  ) +
  scale_x_continuous(
    limits = c(0, 1),
    labels = percent
  ) +
  scale_shape_manual(
    values = c("FALSE" = 16, "TRUE" = 1),
    labels = c("At least 10 samples", "Fewer than 10 samples"),
    name = NULL
  ) +
  labs(
    title = "Assay coverage within sample subgroups",
    x = "Samples represented",
    y = NULL,
    caption = paste(
      "Dashed line marks overall layer coverage.",
      selection_shift_statement
    )
  ) +
  theme(
    strip.text.y = element_text(angle = 0),
    strip.text = element_text(size = 10),
    axis.text = element_text(size = 10),
    panel.grid.major.y = element_blank()
  )
Treatment-response, microarray, and CNV coverage within sex and derived race groups.

Figure 5: Treatment-response, microarray, and CNV coverage within sex and derived race groups.

The full collection is compositionally imbalanced by sex and derived race. Largest subgroup shift among groups with at least 10 samples: Missing in CNV (-45.8% versus overall coverage). Small subgroups have unstable coverage percentages. These are descriptive representation differences; no external reference population is used.

Primary site provides a second view of collection and assay-selection concentration.

bias_samples[,
  primary_site := fcase(
    is.na(ccle_primary_site) |
      trimws(as.character(ccle_primary_site)) == "",
    "Missing",
    default = as.character(ccle_primary_site)
  )
]
top_primary_sites <- bias_samples[
  primary_site != "Missing",
  .N,
  by = primary_site
][
  order(-N, primary_site)
][
  seq_len(min(10L, .N)),
  primary_site
]
bias_samples[,
  primary_site_group := fcase(
    primary_site == "Missing",
    "Missing",
    primary_site %in% top_primary_sites,
    primary_site,
    default = "Other"
  )
]
bias_samples[, metadata := TRUE]
site_layer_membership <- melt(
  bias_samples[, .(
    primary_site_group,
    metadata,
    treatment_response,
    microarray,
    cnv
  )],
  id.vars = "primary_site_group",
  measure.vars = c(
    "metadata",
    "treatment_response",
    "microarray",
    "cnv"
  ),
  variable.name = "data_layer",
  value.name = "present"
)[present == TRUE]
primary_site_composition <- site_layer_membership[
  ,
  .(samples = .N),
  by = .(data_layer, primary_site_group)
]
primary_site_composition[,
  fraction := samples / sum(samples),
  by = data_layer
]
primary_site_composition[,
  data_layer := factor(
    data_layer,
    levels = c("metadata", "treatment_response", "microarray", "cnv"),
    labels = c("Metadata", "Treatment response", "Microarray", "CNV")
  )
]
primary_site_composition[,
  primary_site_group_label := factor(
    gsub("_", " ", primary_site_group, fixed = TRUE),
    levels = gsub(
      "_",
      " ",
      c(top_primary_sites, "Other", "Missing"),
      fixed = TRUE
    )
  )
]
primary_site_colours <- c(
  "#4E79A7",
  "#F28E2B",
  "#E15759",
  "#76B7B2",
  "#59A14F",
  "#EDC948",
  "#B07AA1",
  "#FF9DA7",
  "#9C755F",
  "#17BECF"
)
primary_site_palette <- setNames(
  c(
    primary_site_colours[seq_along(top_primary_sites)],
    "#6B7280",
    "#D1D5DB"
  ),
  gsub(
    "_",
    " ",
    c(top_primary_sites, "Other", "Missing"),
    fixed = TRUE
  )
)

ggplot(
  primary_site_composition,
  aes(x = data_layer, y = fraction, fill = primary_site_group_label)
) +
  geom_col(width = 0.72) +
  scale_y_continuous(labels = percent) +
  scale_fill_manual(
    values = primary_site_palette,
    drop = FALSE,
    name = "Primary site"
  ) +
  labs(
    title = "Primary-site composition across data layers",
    x = NULL,
    y = "Percentage",
    caption = "Ten most frequent metadata sites shown; remaining sites grouped as Other."
  ) +
  guides(fill = guide_legend(ncol = 3))
Primary-site composition of metadata, treatment-response, microarray, and CNV layers.

Figure 6: Primary-site composition of metadata, treatment-response, microarray, and CNV layers.

Lung and haematopoietic/lymphoid samples are the two largest primary-site groups. Site composition reflects the assembled cell-line collection, not cancer incidence.

Molecular Profiles

The loaded PSet contains 2 molecular profiles: microarray.gene_expression_log2 and cnv.gene_log2. No additional molecular profiles are present.

Microarray Expression

The microarray.gene_expression_log2 profile contains 17,970 genes across 823 samples. Its internal assay is exprs. There are 0 non-Ensembl row names, 0 duplicated row names, and 573 source features excluded during mapping. Missing values account for 0.000% of the matrix.

Values are average log2-transformed Affymetrix microarray expression from the CTRPv2.1 source release.

The expression density uses every finite value in the matrix and displays the axis from zero.

microarray_values <- SummarizedExperiment::assay(
  microarray_profile,
  "exprs"
)
microarray_distribution <- data.table(
  expression = as.numeric(microarray_values)
)[is.finite(expression)]
microarray_density <- density(
  microarray_distribution$expression,
  n = 2048,
  from = 0,
  to = max(microarray_distribution$expression)
)
microarray_density <- data.table(
  expression = microarray_density$x,
  density = microarray_density$y
)

ggplot(microarray_density, aes(x = expression, y = density)) +
  geom_area(
    fill = "#4C78A8",
    colour = "#2D4D69",
    alpha = 0.45,
    linewidth = 0.7
  ) +
  geom_vline(
    xintercept = median(microarray_distribution$expression),
    colour = "#9B3A42",
    linetype = "dashed"
  ) +
  scale_x_continuous(
    limits = c(0, max(microarray_distribution$expression)),
    expand = expansion(mult = c(0, 0.02))
  ) +
  scale_y_continuous(expand = expansion(mult = c(0, 0.05))) +
  labs(
    title = "Microarray expression distribution",
    subtitle = paste(
      "n =",
      comma(nrow(microarray_distribution))
    ),
    x = "Average log2-transformed Affymetrix expression",
    y = "Probability density (area = 1)",
    caption = "Dashed line marks median."
  )
Distribution of all finite microarray expression values.

Figure 7: Distribution of all finite microarray expression values.

Copy-Number Variation

The cnv.gene_log2 profile contains 22,598 genes across 823 samples. Its internal assay is copy_number. There are 0 non-Ensembl row names, 0 duplicated row names, and 576 source features excluded during mapping. Missing values account for 0.000% of the matrix.

Most gene-level CNV values should cluster around a neutral log2 ratio of zero, with broader tails representing copy-number losses and gains. The density uses every finite value and a symmetric x-axis centered on zero. Y-axis shows probability density, not copy-number units; area under curve equals 1, so peak height can exceed 1.

cnv_values <- SummarizedExperiment::assay(cnv_profile, "copy_number")
cnv_distribution <- data.table(
  copy_number = as.numeric(cnv_values)
)[is.finite(copy_number)]
cnv_display_limit <- max(abs(range(cnv_distribution$copy_number)))
cnv_density <- density(
  cnv_distribution$copy_number,
  n = 4096,
  from = -cnv_display_limit,
  to = cnv_display_limit
)
cnv_density <- data.table(
  copy_number = cnv_density$x,
  density = cnv_density$y
)

ggplot(cnv_density, aes(x = copy_number, y = density)) +
  geom_area(
    fill = "#59A14F",
    colour = "#356A31",
    alpha = 0.45,
    linewidth = 0.7
  ) +
  geom_vline(
    xintercept = 0,
    colour = "#9B3A42",
    linetype = "dashed"
  ) +
  scale_x_continuous(
    limits = c(-cnv_display_limit, cnv_display_limit),
    expand = expansion(mult = c(0.01, 0.01))
  ) +
  scale_y_continuous(expand = expansion(mult = c(0, 0.05))) +
  labs(
    title = "Copy-number distribution",
    subtitle = paste(
      "n =",
      comma(nrow(cnv_distribution))
    ),
    x = "Gene-level copy-number log2 ratio",
    y = "Probability density (area = 1)",
    caption = "Dashed line marks log2 ratio of zero."
  )
Distribution of all finite gene-level copy-number values.

Figure 8: Distribution of all finite gene-level copy-number values.

Treatment-Response Coverage

CTRPv2 tests each treatment across many cell lines, while the number of treatments observed for a given cell line varies more widely.

All AUC, AAC, relative IC50, absolute AC50, fit-status, and fitted-curve values labeled as recomputed in this report were calculated from the stored CTRPv2 dose-response measurements using PharmacoGx. Recomputed AUC and AAC integrate the PharmacoGx fitted curve. Published CTRP metrics are kept in separate fields and are identified explicitly.

Metric Provenance

Treatment-response profiles stores CTRP-reported values beside PharmacoGx recomputed values:

  • Reported by CTRP: published_auc and published_ec50_umol_raw.
  • Cleaned CTRP comparison field: published_ec50_umol retains positive source EC50 values at or below 1,000,000 µM. This is a cleaned source metric, not a PharmacoGx recomputation.
  • CTRP dose-response series in sensitivity: dose_umol_*, viability_*, predicted_viability_*, and viability_errorbar_*. predicted_viability_* stores CTRP’s source curve predictions at the tested concentrations.
  • CTRP curve-fit fields retained as source provenance: p1_center is log2 apparent EC50, p2_slope is the slope parameter, p4_baseline is the high-concentration asymptote, and p3_total_decline is derived as 1 - p4_baseline. pred_pv_high_conc, conc_pts_fit, and fit_num_param retain the remaining source fit information.
  • PharmacoGx recomputed fit parameters: HS, E0, E_inf, EC50, and Rsq.
  • PharmacoGx recomputed response metrics: auc_recomputed, aac_recomputed, ic50_recomputed, and ac50_recomputed. AUC and AAC use area.type = "Fitted"; both integrate the PharmacoGx curve between the lowest and highest tested concentrations and normalize by that log-dose span.
  • PharmacoGx recomputation QC: recomputed_fit_status, recomputed_fit_message, recomputed_degenerate_fit, recomputed_fitted_crosses_absolute_50, recomputed_empirical_crosses_absolute_50, recomputed_valid_relative_ic50, recomputed_valid_ac50, recomputed_relative_ic50_below_tested_range, recomputed_relative_ic50_above_tested_range, recomputed_ac50_below_tested_range, and recomputed_ac50_above_tested_range.
sample_response_counts <- profiles[,
  .(response_profiles = .N),
  by = master_ccl_id
][, entity := "Samples"]
treatment_response_counts <- profiles[,
  .(response_profiles = .N),
  by = master_cpd_id
][, entity := "Treatments"]
response_counts <- rbindlist(
  list(
    sample_response_counts[, .(response_profiles, entity)],
    treatment_response_counts[, .(response_profiles, entity)]
  )
)

ggplot(response_counts, aes(x = response_profiles)) +
  geom_histogram(
    bins = 35,
    fill = "#4C78A8",
    colour = "white"
  ) +
  facet_wrap(~entity, scales = "free") +
  scale_x_continuous(labels = comma) +
  scale_y_continuous(
    trans = scales::transform_log1p(),
    breaks = c(0, 1, 3, 10, 30, 100, 300, 1000),
    labels = comma,
    minor_breaks = NULL
  ) +
  labs(
    title = "Response-profile coverage",
    x = "Profiles per sample or treatment",
    y = "Records (log[count + 1] scale)",
    caption = paste(
      comma(uniqueN(profiles$master_ccl_id)),
      "of",
      comma(nrow(samples)),
      "samples and",
      comma(uniqueN(profiles$master_cpd_id)),
      "of",
      comma(nrow(treatments)),
      "treatments have response profiles."
    )
  )
Treatment-response profile coverage across samples and treatments.

Figure 9: Treatment-response profile coverage across samples and treatments.

The PSet contains 545 treatment metadata records with 545 unique treatment identifiers. All 545 of those treatments occur in the treatment-response profiles.

The number of finite concentration points indicates how much observed data supported each PharmacoGx curve fit.

fit_point_counts <- profiles[
  is.finite(n_fit_points),
  .(n_fit_points = as.integer(n_fit_points))
]

ggplot(fit_point_counts, aes(x = n_fit_points)) +
  geom_bar(fill = "#72A0C1", width = 0.8) +
  scale_y_log10(labels = comma) +
  labs(
    title = "Observed concentration points per response profile",
    x = "Finite concentration points used for fitting",
    y = "Response profiles",
    caption = paste(
      comma(nrow(fit_point_counts)),
      "profiles have a recorded fit-point count."
    )
  )
Number of observed concentration points supporting each recomputed fit.

Figure 10: Number of observed concentration points supporting each recomputed fit.

PharmacoGx reports both area under the fitted viability curve (AUC) and area above it (AAC). Both use the tested log-concentration range and are complementary views of the same response profiles.

area_metrics <- rbindlist(list(
  profiles[
    is.finite(auc_recomputed),
    .(metric = "Fitted-curve AUC", value = as.numeric(auc_recomputed))
  ],
  profiles[
    is.finite(aac_recomputed),
    .(metric = "Fitted-curve AAC", value = as.numeric(aac_recomputed))
  ]
))

ggplot(area_metrics, aes(x = value)) +
  geom_histogram(
    bins = 60,
    fill = "#6B79B8",
    colour = "white"
  ) +
  facet_wrap(~metric, scales = "free_y") +
  scale_y_log10(labels = comma) +
  labs(
    title = "PharmacoGx fitted-curve area metrics",
    x = "Normalized area (%)",
    y = "Response profiles"
  )
Distributions of PharmacoGx recomputed fitted-curve AUC and AAC.

Figure 11: Distributions of PharmacoGx recomputed fitted-curve AUC and AAC.

Relative IC50 is defined at half of the fitted response span. AC50 is more restrictive: it is defined only when the fitted response crosses absolute 50% viability.

potency_metrics <- rbindlist(list(
  profiles[
    is.finite(ic50_recomputed) & ic50_recomputed > 0,
    .(
      metric = "Relative IC50",
      log10_concentration = log10(as.numeric(ic50_recomputed))
    )
  ],
  profiles[
    is.finite(ac50_recomputed) & ac50_recomputed > 0,
    .(
      metric = "Absolute AC50",
      log10_concentration = log10(as.numeric(ac50_recomputed))
    )
  ]
))

ggplot(potency_metrics, aes(x = log10_concentration)) +
  geom_histogram(
    bins = 60,
    fill = "#D08A4B",
    colour = "white"
  ) +
  facet_wrap(~metric, scales = "free_y") +
  scale_y_log10(labels = comma) +
  labs(
    title = "Recomputed potency metrics",
    x = expression(log[10] * " concentration (" * mu * "M)"),
    y = "Response profiles",
    caption = paste(
      comma(sum(
        is.finite(profiles$ic50_recomputed) &
          profiles$ic50_recomputed > 0
      )),
      "profiles have relative IC50;",
      comma(sum(
        is.finite(profiles$ac50_recomputed) &
          profiles$ac50_recomputed > 0
      )),
      "have absolute AC50."
    )
  )
Distributions of PharmacoGx recomputed relative IC50 and absolute AC50.

Figure 12: Distributions of PharmacoGx recomputed relative IC50 and absolute AC50.

Example Dose-Response Fits

Each plot compares four cell lines treated with the same drug. Points are CTRP observations stored in the PSet. Lines reconstruct the PharmacoGx four-parameter logistic fit from HS, E0, E_inf, and EC50. Highest-R² successful profile is shown when a drug-cell-line pair has replicate profiles.

sensitivity <- as.data.table(CoreGx::assay(treatment_response, "sensitivity"))
example_specification <- data.table(
  Compound.Name = c(
    rep("bortezomib", 4),
    rep("erlotinib", 4)
  ),
  CELL_NAME = c(
    "A375",
    "A549",
    "HCT116",
    "MCF7",
    "HCC827",
    "CALU3",
    "A549",
    "NCIH1975"
  ),
  curve_order = rep(1:4, 2)
)
example_profile_candidates <- merge(
  profiles,
  samples[, .(
    master_ccl_id,
    CELL_NAME = as.character(CELL_NAME)
  )],
  by = "master_ccl_id",
  all.x = TRUE
)
example_profile_candidates <- merge(
  example_profile_candidates,
  example_specification,
  by = c("Compound.Name", "CELL_NAME"),
  all = FALSE
)
example_profiles <- example_profile_candidates[
  recomputed_fit_status == "fit_success" &
    is.finite(HS) &
    is.finite(E0) &
    is.finite(E_inf) &
    is.finite(EC50) &
    EC50 > 0
][
  order(Compound.Name, curve_order, -Rsq, -n_fit_points, response_id),
  .SD[1],
  by = .(Compound.Name, CELL_NAME, curve_order)
]

example_observed <- melt(
  sensitivity[response_id %in% example_profiles$response_id],
  id.vars = c("response_id", "master_ccl_id"),
  measure.vars = patterns(
    dose_umol = "^dose_umol_[0-9]+$",
    viability = "^viability_[0-9]+$"
  ),
  variable.name = "dose_index"
)
example_observed <- merge(
  example_observed,
  example_profiles[, .(
    response_id,
    Compound.Name,
    CELL_NAME,
    curve_order
  )],
  by = "response_id",
  all.x = TRUE
)[
  is.finite(dose_umol) & dose_umol > 0 & is.finite(viability)
]

example_fit_curves <- example_profiles[,
  {
    dose_sequence <- 10^seq(
      log10(dose_min_umol),
      log10(dose_max_umol),
      length.out = 250
    )
    .(
      dose_umol = dose_sequence,
      fitted_viability = E_inf +
        (E0 - E_inf) / (1 + (dose_sequence / EC50)^HS)
    )
  },
  by = .(Compound.Name, CELL_NAME, curve_order)
]
bortezomib_profiles <- example_profiles[
  Compound.Name == "bortezomib"
][order(curve_order)]
bortezomib_cell_lines <- as.character(bortezomib_profiles$CELL_NAME)
bortezomib_observed <- example_observed[
  Compound.Name == "bortezomib"
]
bortezomib_observed[,
  CELL_NAME := factor(
    CELL_NAME,
    levels = bortezomib_cell_lines
  )
]
bortezomib_fits <- example_fit_curves[
  Compound.Name == "bortezomib"
]
bortezomib_fits[,
  CELL_NAME := factor(
    CELL_NAME,
    levels = bortezomib_cell_lines
  )
]
bortezomib_colours <- setNames(
  c("#3B6FB6", "#D17A22", "#3D8B5C", "#9A5BA5"),
  bortezomib_cell_lines
)

ggplot() +
  geom_hline(yintercept = 50, colour = "grey60", linetype = "dotted") +
  geom_line(
    data = bortezomib_fits,
    aes(
      x = dose_umol,
      y = fitted_viability,
      colour = CELL_NAME
    ),
    linewidth = 1.1
  ) +
  geom_point(
    data = bortezomib_observed,
    aes(x = dose_umol, y = viability, colour = CELL_NAME),
    fill = "white",
    shape = 21,
    size = 2.4,
    stroke = 0.8
  ) +
  scale_colour_manual(values = bortezomib_colours) +
  scale_x_log10(labels = label_number()) +
  coord_cartesian(ylim = c(0, 120)) +
  labs(
    title = "Bortezomib response across four cell lines",
    subtitle = "Points: CTRP observations; lines: PharmacoGx recomputed fits",
    x = expression("Concentration (" * mu * "M, log scale)"),
    y = "Viability (%)",
    colour = "Cell line",
    caption = "Dotted line marks absolute 50% viability."
  )
Bortezomib dose-response observations and PharmacoGx recomputed fits for four cell lines.

Figure 13: Bortezomib dose-response observations and PharmacoGx recomputed fits for four cell lines. 

erlotinib_profiles <- example_profiles[
  Compound.Name == "erlotinib"
][order(curve_order)]
erlotinib_cell_lines <- as.character(erlotinib_profiles$CELL_NAME)
erlotinib_observed <- example_observed[
  Compound.Name == "erlotinib"
]
erlotinib_observed[,
  CELL_NAME := factor(
    CELL_NAME,
    levels = erlotinib_cell_lines
  )
]
erlotinib_fits <- example_fit_curves[
  Compound.Name == "erlotinib"
]
erlotinib_fits[,
  CELL_NAME := factor(
    CELL_NAME,
    levels = erlotinib_cell_lines
  )
]
erlotinib_colours <- setNames(
  c("#3B6FB6", "#D17A22", "#3D8B5C", "#9A5BA5"),
  erlotinib_cell_lines
)

ggplot() +
  geom_hline(yintercept = 50, colour = "grey60", linetype = "dotted") +
  geom_line(
    data = erlotinib_fits,
    aes(
      x = dose_umol,
      y = fitted_viability,
      colour = CELL_NAME
    ),
    linewidth = 1.1
  ) +
  geom_point(
    data = erlotinib_observed,
    aes(x = dose_umol, y = viability, colour = CELL_NAME),
    fill = "white",
    shape = 21,
    size = 2.4,
    stroke = 0.8
  ) +
  scale_colour_manual(values = erlotinib_colours) +
  scale_x_log10(labels = label_number()) +
  coord_cartesian(ylim = c(0, 120)) +
  labs(
    title = "Erlotinib response across four cell lines",
    subtitle = "Points: CTRP observations; lines: PharmacoGx recomputed fits",
    x = expression("Concentration (" * mu * "M, log scale)"),
    y = "Viability (%)",
    colour = "Cell line",
    caption = "Dotted line marks absolute 50% viability."
  )
Erlotinib dose-response observations and PharmacoGx recomputed fits for four cell lines.

Figure 14: Erlotinib dose-response observations and PharmacoGx recomputed fits for four cell lines.

Published CTRP vs PharmacoGx Metrics

Rees et al. fit CTRP’s 16-point concentration-response series with a three-parameter sigmoid whose low-concentration asymptote is fixed at 100% viability. When the initial apparent EC50 exceeded the highest tested concentration, CTRP refit the curve with a two-parameter sigmoid whose high-concentration asymptote is fixed at 0%. The PSet contains 306,346 three-parameter and 88,917 two-parameter source fits.

Published apparent EC50 is the concentration at half of the CTRP fitted response span. In a two-parameter source fit it is also an absolute 50% viability crossing because the source asymptotes are fixed at 100% and 0%. PharmacoGx relative IC50 is the midpoint of its independently fitted response span. PharmacoGx AC50 instead represents a crossing of absolute 50% viability.

Raw published EC50 values are retained in the PSet for provenance. The cleaned field removes nonpositive values and values above the configured 1,000,000 µM cutoff before comparative analysis.

published_ec50_cutoff_umol <- 1000000
ec50_retention <- profiles[,
  .(
    category = fcase(
      is.finite(published_ec50_umol)                                                            ,
      "Retained for comparison"                                                                 ,
      is.finite(published_ec50_umol_raw) & published_ec50_umol_raw <= 0                         ,
      "Removed: nonpositive"                                                                    ,
      is.finite(published_ec50_umol_raw) & published_ec50_umol_raw > published_ec50_cutoff_umol ,
      paste0("Removed: above ", comma(published_ec50_cutoff_umol), " µM cutoff")                ,
      default = "Missing raw value"
    )
  )
][,
  .(response_profiles = .N),
  by = category
][order(response_profiles)]

ggplot(
  ec50_retention,
  aes(x = response_profiles, y = reorder(category, response_profiles))
) +
  geom_col(fill = "#B26B3F", width = 0.7) +
  geom_text(
    aes(label = comma(response_profiles)),
    hjust = -0.12,
    colour = "grey20"
  ) +
  scale_x_continuous(
    labels = comma,
    expand = expansion(mult = c(0, 0.14))
  ) +
  labs(
    title = "Published EC50 values available for comparison",
    x = "Response profiles",
    y = NULL
  )
Retention of published CTRP apparent EC50 values for comparison.

Figure 15: Retention of published CTRP apparent EC50 values for comparison.

Published CTRP AUC integrates its fitted sigmoid over the standard 16-point, twofold concentration series. This design contains 15 log2 concentration intervals, so a source AUC near 15 corresponds approximately to 100% viability across the assay range.

PharmacoGx auc_recomputed uses computeAUC() with area.type = "Fitted". It integrates the PharmacoGx four-parameter curve over the tested log-concentration range, divides by that range, and reports normalized average fitted viability on a percentage scale. aac_recomputed is the complementary area above the same fitted curve.

For comparison, published AUC is converted to the PharmacoGx scale with published_auc * 100 / 15. This conversion aligns units and makes the comparison fit-to-fit, but it does not make the estimators identical. CTRP uses its constrained three-parameter or conditional two-parameter model over the standard assay design. PharmacoGx uses an independently fitted four-parameter curve over each profile’s available tested range. CTRP does not provide a separate normalization formula, so the factor of 15 is derived from the documented 16-point, twofold design.

auc_pairs <- profiles[
  is.finite(published_auc) & is.finite(auc_recomputed),
  .(
    published_auc = as.numeric(published_auc),
    published_auc_standardized = as.numeric(published_auc) * 100 / 15,
    recomputed_auc = as.numeric(auc_recomputed),
    fit_num_param = as.integer(fit_num_param)
  )
]
auc_pairs[,
  source_fit := factor(
    paste0("CTRP ", fit_num_param, "-parameter fit"),
    levels = c("CTRP 2-parameter fit", "CTRP 3-parameter fit")
  )
]
auc_pearson <- cor(
  auc_pairs$published_auc_standardized,
  auc_pairs$recomputed_auc
)
auc_spearman <- cor(
  auc_pairs$published_auc_standardized,
  auc_pairs$recomputed_auc,
  method = "spearman"
)
auc_median_difference <- median(
  auc_pairs$recomputed_auc - auc_pairs$published_auc_standardized
)
auc_median_absolute_difference <- median(
  abs(auc_pairs$recomputed_auc - auc_pairs$published_auc_standardized)
)
auc_published_above_100 <- auc_pairs[
  published_auc_standardized > 100,
  .N
]
auc_published_above_100_three_parameter <- auc_pairs[
  published_auc_standardized > 100 & fit_num_param == 3,
  .N
]

ggplot(
  auc_pairs,
  aes(x = published_auc_standardized, y = recomputed_auc)
) +
  geom_abline(
    intercept = 0,
    slope = 1,
    colour = "grey45",
    linetype = "dashed",
    linewidth = 0.8
  ) +
  geom_bin_2d(bins = 70) +
  geom_smooth(
    method = "lm",
    formula = y ~ x,
    se = FALSE,
    colour = "#9B3A42",
    linewidth = 0.9
  ) +
  scale_fill_viridis_c(
    trans = "log10",
    labels = comma,
    name = "Profiles per bin"
  ) +
  facet_wrap(~source_fit) +
  coord_cartesian(ylim = range(auc_pairs$recomputed_auc)) +
  labs(
    title = "Published CTRP AUC vs PharmacoGx fitted-curve AUC",
    subtitle = paste0(
      "N = ",
      comma(nrow(auc_pairs)),
      ", Pearson r = ",
      number(auc_pearson, accuracy = 0.001),
      ", Spearman rho = ",
      number(auc_spearman, accuracy = 0.001),
      ", Median absolute difference = ",
      number(auc_median_absolute_difference, accuracy = 0.01),
      " points"
    ),
    x = "Converted published CTRP AUC (%)",
    y = "PharmacoGx fitted-curve AUC (%)",
    caption = paste0(
      "Published AUC × 100 / 15. Dashed line: identity; red line: fitted trend.\n",
      comma(auc_published_above_100),
      " converted CTRP values exceed 100; ",
      comma(auc_published_above_100_three_parameter),
      " come from three-parameter source fits."
    )
  )
Converted published CTRP AUC compared with PharmacoGx recomputed fitted-curve AUC.

Figure 16: Converted published CTRP AUC compared with PharmacoGx recomputed fitted-curve AUC.

After conversion, median recomputed-minus-published difference is -0.98 percentage points. Matched distributions below use the same x-axis.

auc_distributions <- rbindlist(list(
  auc_pairs[, .(
    metric = "Published CTRP AUC (converted)",
    auc = published_auc_standardized
  )],
  auc_pairs[, .(
    metric = "PharmacoGx fitted-curve AUC",
    auc = recomputed_auc
  )]
))
auc_distributions[,
  metric := factor(
    metric,
    levels = c(
      "Published CTRP AUC (converted)",
      "PharmacoGx fitted-curve AUC"
    )
  )
]

ggplot(auc_distributions, aes(x = auc, fill = metric)) +
  geom_histogram(
    bins = 60,
    colour = "white",
    linewidth = 0.25
  ) +
  facet_wrap(~metric, ncol = 1) +
  scale_fill_manual(
    values = c(
      "Published CTRP AUC (converted)" = "#D17A22",
      "PharmacoGx fitted-curve AUC" = "#3B6FB6"
    ),
    guide = "none"
  ) +
  scale_y_continuous(labels = comma) +
  labs(
    title = "Published and PharmacoGx fitted AUC distributions",
    subtitle = paste(comma(nrow(auc_pairs)), "matched response profiles"),
    x = "AUC (normalized viability %)",
    y = "Response profiles",
    caption = paste(
      "Published conversion: AUC × 100 / 15."
    )
  )
Matched published CTRP and PharmacoGx recomputed AUC distributions.

Figure 17: Matched published CTRP and PharmacoGx recomputed AUC distributions.

Both distributions concentrate near high viability with long lower tails. Converted CTRP AUC can exceed 100 when the source three-parameter fit has a negative total decline and a high-concentration asymptote above the 100% control level.

Published CTRP vs PharmacoGx Fit Example

Comparison uses one response profile shared by the source and recomputed data.

  • Observed points: CTRP viability_* values paired with dose_umol_*.
  • CTRP published fit: source-supplied predicted_viability_* values paired with the same numbered doses, ordered by concentration, then connected with a dashed line. No CTRP curve is refitted here, and the line is not a continuous reconstruction from p1_center, p2_slope, p3_total_decline, or p4_baseline. fit_num_param identifies whether CTRP retained its three-parameter fit or used the conditional two-parameter refit.
  • PharmacoGx recomputed fit: stored E0, E_inf, EC50, and HS parameters evaluated with the PharmacoGx four-parameter Hill equation at 300 log-spaced concentrations between the lowest and highest displayed CTRP doses.

Both lines therefore represent the same drug, cell line, response identifier, and dose range. Difference: CTRP line displays source predictions only at the tested doses; PharmacoGx line is a continuous reconstruction of the recomputed fit.

curve_comparison_profile <- merge(
  profiles,
  samples[, .(
    master_ccl_id,
    CELL_NAME = as.character(CELL_NAME)
  )],
  by = "master_ccl_id",
  all.x = TRUE
)[
  Compound.Name == "erlotinib" &
    CELL_NAME == "HCC827" &
    recomputed_fit_status == "fit_success" &
    is.finite(HS) &
    is.finite(E0) &
    is.finite(E_inf) &
    is.finite(EC50) &
    EC50 > 0
][
  order(-Rsq, -n_fit_points, response_id)
][1]

curve_comparison_points <- melt(
  sensitivity[response_id == curve_comparison_profile$response_id],
  id.vars = c("response_id", "master_ccl_id"),
  measure.vars = list(
    dose_umol = grep("^dose_umol_[0-9]+$", names(sensitivity)),
    observed_viability = grep("^viability_[0-9]+$", names(sensitivity)),
    ctrp_predicted_viability = grep(
      "^predicted_viability_[0-9]+$",
      names(sensitivity)
    )
  ),
  variable.name = "dose_index"
)[
  is.finite(dose_umol) &
    dose_umol > 0 &
    is.finite(observed_viability) &
    is.finite(ctrp_predicted_viability)
][order(dose_umol)]

curve_comparison_points[,
  pgx_predicted_viability := curve_comparison_profile$E_inf +
    (curve_comparison_profile$E0 - curve_comparison_profile$E_inf) /
      (1 +
        (dose_umol /
          curve_comparison_profile$EC50)^curve_comparison_profile$HS)
]

curve_comparison_dose_sequence <- 10^seq(
  log10(min(curve_comparison_points$dose_umol)),
  log10(max(curve_comparison_points$dose_umol)),
  length.out = 300
)
curve_comparison_pharmacogx <- data.table(
  dose_umol = curve_comparison_dose_sequence,
  fitted_viability = curve_comparison_profile$E_inf +
    (curve_comparison_profile$E0 - curve_comparison_profile$E_inf) /
      (1 +
        (curve_comparison_dose_sequence /
          curve_comparison_profile$EC50)^curve_comparison_profile$HS),
  fit_source = "PharmacoGx recomputed fit"
)
curve_comparison_ctrp <- curve_comparison_points[, .(
  dose_umol,
  fitted_viability = ctrp_predicted_viability,
  fit_source = "CTRP published fit"
)]
curve_comparison_fits <- rbindlist(list(
  curve_comparison_ctrp,
  curve_comparison_pharmacogx
))
curve_comparison_fits[,
  fit_source := factor(
    fit_source,
    levels = c("CTRP published fit", "PharmacoGx recomputed fit")
  )
]

ggplot() +
  geom_hline(yintercept = 50, colour = "grey65", linetype = "dotted") +
  geom_line(
    data = curve_comparison_fits,
    aes(
      x = dose_umol,
      y = fitted_viability,
      colour = fit_source,
      linetype = fit_source
    ),
    linewidth = 1.15
  ) +
  geom_point(
    data = curve_comparison_points,
    aes(x = dose_umol, y = observed_viability),
    shape = 21,
    size = 2.6,
    stroke = 0.8,
    colour = "grey20",
    fill = "white"
  ) +
  scale_colour_manual(
    values = c(
      "CTRP published fit" = "#D17A22",
      "PharmacoGx recomputed fit" = "#3B6FB6"
    )
  ) +
  scale_linetype_manual(
    values = c(
      "CTRP published fit" = "dashed",
      "PharmacoGx recomputed fit" = "solid"
    )
  ) +
  scale_x_log10(labels = label_number()) +
  coord_cartesian(ylim = c(0, 120)) +
  labs(
    title = "Erlotinib — HCC827 fitted-curve comparison",
    subtitle = paste0(
      "Points: observed CTRP viability; source fit: ",
      curve_comparison_profile$fit_num_param,
      "-parameter"
    ),
    x = expression("Concentration (" * mu * "M, log scale)"),
    y = "Viability (%)",
    colour = NULL,
    linetype = NULL,
    caption = paste0(
      "Response profile: ",
      curve_comparison_profile$response_id,
      "."
    )
  )
Published CTRP predictions and PharmacoGx recomputed fit for one response profile.

Figure 18: Published CTRP predictions and PharmacoGx recomputed fit for one response profile.

Paired Fit-to-Point Comparison

Both fits can be scored against the same observed viability values at the same concentrations. Evaluation R² is calculated as 1 - sum((observed - predicted)^2) / sum((observed - mean(observed))^2). This is calculated fresh for both sources; stored source-specific R² fields are not compared.

RMSE is the primary comparison because it remains in viability percentage points. MAE is included as a less outlier-sensitive error measure. Higher R² is better; lower RMSE and MAE are better. For one response profile, R² and RMSE rank the fits identically because both use the same observed points and differ only through squared error.

curve_comparison_observed_mean <- mean(
  curve_comparison_points$observed_viability
)
curve_comparison_observed_sst <- sum(
  (curve_comparison_points$observed_viability -
    curve_comparison_observed_mean)^2
)

curve_comparison_scores <- data.table(
  fit_source = c("CTRP published fit", "PharmacoGx recomputed fit"),
  r_squared = c(
    1 -
      sum(
        (curve_comparison_points$observed_viability -
          curve_comparison_points$ctrp_predicted_viability)^2
      ) /
        curve_comparison_observed_sst,
    1 -
      sum(
        (curve_comparison_points$observed_viability -
          curve_comparison_points$pgx_predicted_viability)^2
      ) /
        curve_comparison_observed_sst
  ),
  rmse = c(
    sqrt(mean(
      (curve_comparison_points$observed_viability -
        curve_comparison_points$ctrp_predicted_viability)^2
    )),
    sqrt(mean(
      (curve_comparison_points$observed_viability -
        curve_comparison_points$pgx_predicted_viability)^2
    ))
  ),
  mae = c(
    mean(abs(
      curve_comparison_points$observed_viability -
        curve_comparison_points$ctrp_predicted_viability
    )),
    mean(abs(
      curve_comparison_points$observed_viability -
        curve_comparison_points$pgx_predicted_viability
    ))
  )
)

curve_comparison_scores_long <- melt(
  curve_comparison_scores,
  id.vars = "fit_source",
  measure.vars = c("r_squared", "rmse", "mae"),
  variable.name = "metric",
  value.name = "value"
)
curve_comparison_scores_long[,
  metric := factor(
    metric,
    levels = c("r_squared", "rmse", "mae"),
    labels = c(
      "Evaluation R² — higher is better",
      "RMSE — lower is better",
      "MAE — lower is better"
    )
  )
]
curve_comparison_scores_long[,
  label := fifelse(
    metric == "Evaluation R² — higher is better",
    number(value, accuracy = 0.001),
    paste0(number(value, accuracy = 0.01), " points")
  )
]

ggplot(
  curve_comparison_scores_long,
  aes(x = value, y = fit_source, fill = fit_source)
) +
  geom_col(width = 0.62) +
  geom_text(
    aes(label = label),
    hjust = -0.12,
    colour = "grey20",
    size = 3.6
  ) +
  facet_wrap(
    ~metric,
    ncol = 1,
    scales = "free_x"
  ) +
  scale_fill_manual(
    values = c(
      "CTRP published fit" = "#D17A22",
      "PharmacoGx recomputed fit" = "#3B6FB6"
    ),
    guide = "none"
  ) +
  scale_x_continuous(
    expand = expansion(mult = c(0, 0.2))
  ) +
  labs(
    title = "Fit-to-point scores for erlotinib — HCC827",
    subtitle = paste(
      nrow(curve_comparison_points),
      "shared observed dose-response points"
    ),
    x = NULL,
    y = NULL
  )
Paired fit-to-point scores for the example response profile.

Figure 19: Paired fit-to-point scores for the example response profile.

Example PharmacoGx fit has higher evaluation R² and lower pointwise error. This does not establish which source performs better across CTRPv2, so the same paired calculation is repeated for every response profile with at least four shared finite points.

fit_comparison_dose_columns <- grep(
  "^dose_umol_[0-9]+$",
  names(sensitivity),
  value = TRUE
)
fit_comparison_observed_columns <- grep(
  "^viability_[0-9]+$",
  names(sensitivity),
  value = TRUE
)
fit_comparison_ctrp_columns <- grep(
  "^predicted_viability_[0-9]+$",
  names(sensitivity),
  value = TRUE
)

fit_comparison_parameters <- profiles[
  match(sensitivity$response_id, response_id),
  .(HS, E0, E_inf, EC50, fit_num_param)
]
fit_comparison_dose_matrix <- as.matrix(
  sensitivity[, ..fit_comparison_dose_columns]
)
fit_comparison_observed_matrix <- as.matrix(
  sensitivity[, ..fit_comparison_observed_columns]
)
fit_comparison_ctrp_matrix <- as.matrix(
  sensitivity[, ..fit_comparison_ctrp_columns]
)
fit_comparison_pgx_matrix <- fit_comparison_parameters$E_inf +
  (fit_comparison_parameters$E0 -
    fit_comparison_parameters$E_inf) /
    (1 +
      (fit_comparison_dose_matrix /
        fit_comparison_parameters$EC50)^fit_comparison_parameters$HS)

fit_comparison_valid_matrix <-
  is.finite(fit_comparison_dose_matrix) &
  fit_comparison_dose_matrix > 0 &
  is.finite(fit_comparison_observed_matrix) &
  is.finite(fit_comparison_ctrp_matrix) &
  is.finite(fit_comparison_pgx_matrix)
fit_comparison_n_points <- rowSums(fit_comparison_valid_matrix)
fit_comparison_observed_mean <- rowSums(ifelse(
  fit_comparison_valid_matrix,
  fit_comparison_observed_matrix,
  0
)) /
  fit_comparison_n_points
fit_comparison_ctrp_sse <- rowSums(ifelse(
  fit_comparison_valid_matrix,
  (fit_comparison_observed_matrix -
    fit_comparison_ctrp_matrix)^2,
  0
))
fit_comparison_pgx_sse <- rowSums(ifelse(
  fit_comparison_valid_matrix,
  (fit_comparison_observed_matrix -
    fit_comparison_pgx_matrix)^2,
  0
))
fit_comparison_observed_sst <- rowSums(ifelse(
  fit_comparison_valid_matrix,
  (fit_comparison_observed_matrix -
    fit_comparison_observed_mean)^2,
  0
))

fit_comparison_scores <- data.table(
  response_id = sensitivity$response_id,
  fit_num_param = fit_comparison_parameters$fit_num_param,
  n_points = fit_comparison_n_points,
  ctrp_r_squared = 1 -
    fit_comparison_ctrp_sse / fit_comparison_observed_sst,
  pgx_r_squared = 1 -
    fit_comparison_pgx_sse / fit_comparison_observed_sst,
  ctrp_rmse = sqrt(
    fit_comparison_ctrp_sse / fit_comparison_n_points
  ),
  pgx_rmse = sqrt(
    fit_comparison_pgx_sse / fit_comparison_n_points
  )
)[
  n_points >= 4 &
    is.finite(ctrp_r_squared) &
    is.finite(pgx_r_squared) &
    is.finite(ctrp_rmse) &
    is.finite(pgx_rmse)
]
fit_comparison_scores[,
  lower_rmse := fcase(
    pgx_rmse < ctrp_rmse - 1e-8 ,
    "PharmacoGx lower error"    ,
    ctrp_rmse < pgx_rmse - 1e-8 ,
    "CTRP lower error"          ,
    default = "Tie"
  )
]
fit_comparison_scores[,
  lower_rmse := factor(
    lower_rmse,
    levels = c(
      "CTRP lower error",
      "Tie",
      "PharmacoGx lower error"
    )
  )
]
fit_comparison_outcomes <- fit_comparison_scores[,
  .(response_profiles = .N),
  by = lower_rmse
]
fit_comparison_outcomes[,
  fraction := response_profiles / sum(response_profiles)
]
fit_comparison_outcomes[, source_fit := "All source fits"]

fit_comparison_outcomes_by_source <- fit_comparison_scores[,
  .(response_profiles = .N),
  by = .(fit_num_param, lower_rmse)
]
fit_comparison_outcomes_by_source[,
  fraction := response_profiles / sum(response_profiles),
  by = fit_num_param
]
fit_comparison_outcomes_by_source[,
  source_fit := paste0("CTRP ", fit_num_param, "-parameter fits")
]
fit_comparison_outcomes_display <- rbindlist(
  list(
    fit_comparison_outcomes[, .(
      source_fit,
      lower_rmse,
      response_profiles,
      fraction
    )],
    fit_comparison_outcomes_by_source[, .(
      source_fit,
      lower_rmse,
      response_profiles,
      fraction
    )]
  )
)
fit_comparison_outcomes_display[,
  source_fit := factor(
    source_fit,
    levels = c(
      "CTRP 3-parameter fits",
      "CTRP 2-parameter fits",
      "All source fits"
    )
  )
]

fit_comparison_median_r_squared_difference <- median(
  fit_comparison_scores$pgx_r_squared -
    fit_comparison_scores$ctrp_r_squared
)
fit_comparison_median_rmse_difference <- median(
  fit_comparison_scores$pgx_rmse -
    fit_comparison_scores$ctrp_rmse
)

ggplot(
  fit_comparison_outcomes_display,
  aes(x = source_fit, y = fraction, fill = lower_rmse)
) +
  geom_col(width = 0.65) +
  geom_text(
    aes(
      label = ifelse(
        fraction >= 0.03,
        paste0(
          comma(response_profiles),
          "\n",
          percent(fraction, accuracy = 0.1)
        ),
        ""
      )
    ),
    position = position_stack(vjust = 0.5),
    colour = "white",
    fontface = "bold",
    lineheight = 1
  ) +
  coord_flip() +
  scale_y_continuous(
    labels = percent,
    limits = c(0, 1),
    expand = expansion(mult = c(0.01, 0.01))
  ) +
  scale_fill_manual(
    values = c(
      "CTRP lower error" = "#D17A22",
      "PharmacoGx lower error" = "#3B6FB6",
      "Tie" = "#8A8F98"
    ),
    breaks = c(
      "PharmacoGx lower error",
      "Tie",
      "CTRP lower error"
    )
  ) +
  labs(
    title = "Which fitted curve is closer to the empirical observations?",
    subtitle = "Paired RMSE using identical points; grouped by CTRP source model",
    x = NULL,
    y = "Fraction of matched response profiles",
    fill = NULL,
    caption = paste0(
      "Tie: ",
      comma(fit_comparison_outcomes[lower_rmse == "Tie", response_profiles]),
      " (",
      percent(
        fit_comparison_outcomes[lower_rmse == "Tie", fraction],
        accuracy = 0.1
      ),
      "). RMSE differences within 1e-8 are treated as ties."
    )
  )
Paired fit-to-point score comparison across eligible response profiles.

Figure 20: Paired fit-to-point score comparison across eligible response profiles.

Across 395,263 matched profiles, CTRP has lower RMSE for 61.4% of matched profiles; PharmacoGx has lower RMSE for 37.9%. Median paired PharmacoGx-minus-CTRP R² difference is -0.002; median paired PharmacoGx-minus-CTRP RMSE difference is 0.10 viability points. CTRP has lower pointwise error more often, but the typical difference is small.

These are in-sample fit-to-point scores. Both curves were estimated from these measurements, so this comparison measures adherence to the available points, not performance at unseen concentrations. CTRP uses a two- or three-parameter source model, while PharmacoGx uses a four-parameter model; the score does not penalize that difference in complexity. Held-out doses or independent replicate experiments would be needed to compare predictive performance.

The potency comparison uses positive cleaned published values and positive recomputed relative IC50 values. Both axes use micromolar concentration on a log10 scale.

ec50_ic50_pairs <- profiles[
  is.finite(published_ec50_umol) &
    published_ec50_umol > 0 &
    is.finite(ic50_recomputed) &
    ic50_recomputed > 0,
  .(
    published_ec50_umol = as.numeric(published_ec50_umol),
    ic50_recomputed = as.numeric(ic50_recomputed),
    log10_published_ec50 = log10(as.numeric(published_ec50_umol)),
    log10_ic50_recomputed = log10(as.numeric(ic50_recomputed))
  )
]
ec50_ic50_log_pearson <- cor(
  ec50_ic50_pairs$log10_published_ec50,
  ec50_ic50_pairs$log10_ic50_recomputed
)
ec50_ic50_spearman <- cor(
  ec50_ic50_pairs$published_ec50_umol,
  ec50_ic50_pairs$ic50_recomputed,
  method = "spearman"
)
ec50_ic50_median_log_difference <- median(abs(
  ec50_ic50_pairs$log10_ic50_recomputed -
    ec50_ic50_pairs$log10_published_ec50
))
ec50_display_x_limits <- quantile(
  ec50_ic50_pairs$log10_published_ec50,
  probs = c(0.05, 0.99)
)
ec50_display_y_limits <- quantile(
  ec50_ic50_pairs$log10_ic50_recomputed,
  probs = c(0.01, 0.99)
)

ggplot(
  ec50_ic50_pairs,
  aes(x = log10_published_ec50, y = log10_ic50_recomputed)
) +
  geom_bin_2d(binwidth = c(0.15, 0.15)) +
  geom_abline(
    slope = 1,
    intercept = 0,
    colour = "#9B3A42",
    linewidth = 0.9
  ) +
  scale_fill_viridis_c(
    trans = "log10",
    labels = comma,
    name = "Profiles per bin"
  ) +
  coord_cartesian(
    xlim = ec50_display_x_limits,
    ylim = ec50_display_y_limits
  ) +
  labs(
    title = "Published apparent EC50 vs recomputed relative IC50",
    subtitle = paste0(
      "N = ",
      comma(nrow(ec50_ic50_pairs)),
      "  |  log10 Pearson r = ",
      number(ec50_ic50_log_pearson, accuracy = 0.001),
      "  |  Spearman rho = ",
      number(ec50_ic50_spearman, accuracy = 0.001),
      "  |  median absolute log10 difference = ",
      number(ec50_ic50_median_log_difference, accuracy = 0.001)
    ),
    x = expression(log[10] * " published apparent EC50 (" * mu * "M)"),
    y = expression(log[10] * " recomputed relative IC50 (" * mu * "M)"),
    caption = paste0(
      "Binning is performed directly in log10 concentration space; the red ",
      "line represents numerical identity.\n",
      "The display shows the central published EC50 and recomputed IC50 ",
      "ranges. Statistics use all positive cleaned pairs."
    )
  )
Published apparent EC50 compared with PharmacoGx recomputed relative IC50.

Figure 21: Published apparent EC50 compared with PharmacoGx recomputed relative IC50.

Recomputed Fit Quality

Fit status, endpoint availability, and range flags answer different questions. A degenerate fit has effectively no fitted response span; a valid relative IC50 can still lie outside the tested concentrations; and AC50 additionally requires the fitted curve to cross absolute 50% viability.

Flags used below mean:

  • Fit success: the curve optimizer returned a nonzero fitted response span.
  • Degenerate fit: the fitted upper and lower responses are effectively the same, so relative IC50 and AC50 are undefined.
  • Too few points: fewer than the configured minimum number of finite, distinct concentration points were available; no fit is reported.
  • Fit failed: the optimizer did not return a usable curve.
  • Valid relative IC50: a finite positive concentration was found at half of the fitted response span.
  • Valid absolute AC50: a finite positive concentration was found where the fitted curve crosses absolute 50% viability.
  • Below/above tested range: the endpoint exists mathematically but lies outside the observed concentration range and is therefore extrapolated.
  • Empirical crossing: the observed viability values bracket 50%; fitted crossing means the fitted curve’s upper and lower response levels bracket 50%.

The current PSet contains only successful and degenerate fit statuses; the other status definitions describe outcomes that the recomputation pipeline can emit when fitting is not possible.

fit_outcomes <- data.table(
  outcome = factor(
    c(
      "Fit success",
      "Degenerate fit",
      "Valid relative IC50",
      "Relative IC50 below tested range",
      "Relative IC50 above tested range",
      "Valid absolute AC50",
      "AC50 below tested range",
      "AC50 above tested range"
    ),
    levels = rev(c(
      "Fit success",
      "Degenerate fit",
      "Valid relative IC50",
      "Relative IC50 below tested range",
      "Relative IC50 above tested range",
      "Valid absolute AC50",
      "AC50 below tested range",
      "AC50 above tested range"
    ))
  ),
  group = c(
    "Fit status",
    "Fit status",
    "Endpoint availability",
    "Range flag",
    "Range flag",
    "Endpoint availability",
    "Range flag",
    "Range flag"
  ),
  response_profiles = c(
    sum(profiles$recomputed_fit_status == "fit_success", na.rm = TRUE),
    sum(profiles$recomputed_degenerate_fit %in% TRUE, na.rm = TRUE),
    sum(profiles$recomputed_valid_relative_ic50 %in% TRUE, na.rm = TRUE),
    sum(
      profiles$recomputed_relative_ic50_below_tested_range %in% TRUE,
      na.rm = TRUE
    ),
    sum(
      profiles$recomputed_relative_ic50_above_tested_range %in% TRUE,
      na.rm = TRUE
    ),
    sum(profiles$recomputed_valid_ac50 %in% TRUE, na.rm = TRUE),
    sum(profiles$recomputed_ac50_below_tested_range %in% TRUE, na.rm = TRUE),
    sum(profiles$recomputed_ac50_above_tested_range %in% TRUE, na.rm = TRUE)
  )
)
fit_outcomes[, fraction := response_profiles / nrow(profiles)]

ggplot(
  fit_outcomes,
  aes(
    x = fraction,
    y = outcome,
    fill = group
  )
) +
  geom_col(width = 0.7) +
  geom_text(
    aes(
      label = paste0(
        comma(response_profiles),
        " (",
        percent(fraction, accuracy = 0.1),
        ")"
      )
    ),
    hjust = -0.08,
    colour = "grey20",
    size = 3.4
  ) +
  scale_x_continuous(
    labels = percent,
    breaks = c(0, 0.25, 0.5, 0.75, 1),
    limits = c(0, 1.15),
    expand = expansion(mult = c(0, 0))
  ) +
  scale_fill_manual(
    values = c(
      "Fit status" = "#4C78A8",
      "Endpoint availability" = "#59A14F",
      "Range flag" = "#E09C45"
    )
  ) +
  labs(
    title = "Fit and endpoint outcomes",
    subtitle = paste(comma(nrow(profiles)), "response profiles"),
    x = "Fraction of response profiles",
    y = NULL,
    fill = NULL
  )
Outcome frequencies for recomputed fits and potency endpoints.

Figure 22: Outcome frequencies for recomputed fits and potency endpoints.

Fitted and empirical crossing calls are shown separately because a fitted curve can cross 50% viability through interpolation or extrapolation even when the observed concentrations do not straddle 50%.

absolute_50_crossing <- profiles[,
  .(response_profiles = .N),
  by = .(
    fitted = factor(
      ifelse(
        recomputed_fitted_crosses_absolute_50 %in% TRUE,
        "Crosses 50%",
        "Does not cross 50%"
      ),
      levels = c("Does not cross 50%", "Crosses 50%")
    ),
    empirical = factor(
      ifelse(
        recomputed_empirical_crosses_absolute_50 %in% TRUE,
        "Crosses 50%",
        "Does not cross 50%"
      ),
      levels = c("Does not cross 50%", "Crosses 50%")
    )
  )
]
absolute_50_crossing[, fraction := response_profiles / nrow(profiles)]

ggplot(
  absolute_50_crossing,
  aes(x = empirical, y = fitted, fill = fraction)
) +
  geom_tile(colour = "white", linewidth = 1.2) +
  geom_text(
    aes(
      label = paste0(
        comma(response_profiles),
        "\n",
        percent(fraction, accuracy = 0.1)
      ),
      colour = fraction > 0.25
    ),
    fontface = "bold",
    size = 4
  ) +
  scale_fill_gradient(
    low = "#E7F0F5",
    high = "#2D6584",
    labels = percent,
    name = "Fraction of\nall profiles"
  ) +
  scale_colour_manual(
    values = c("FALSE" = "grey10", "TRUE" = "white"),
    guide = "none"
  ) +
  labs(
    title = "Fitted versus empirical absolute-50 crossing",
    x = "Observed dose-response points",
    y = "Recomputed fitted curve"
  ) +
  theme(
    legend.position = "right",
    panel.grid = element_blank()
  )
Agreement between fitted and empirical absolute-50 crossing calls.

Figure 23: Agreement between fitted and empirical absolute-50 crossing calls.

Session Info

sessionInfo()
## R version 4.4.3 (2025-02-28)
## Platform: aarch64-apple-darwin20.0.0
## Running under: macOS 27.0
## 
## Matrix products: default
## BLAS/LAPACK: /Users/michael/Projects/BHKLab/ARPA-H/CTRP-Pharmacoset_Snakemake_new/.pixi/envs/default/lib/libopenblas.0.dylib;  LAPACK version 3.12.0
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## time zone: America/Toronto
## tzcode source: system (macOS)
## 
## attached base packages:
## [1] stats4    stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] MultiAssayExperiment_1.32.0 PharmacoGx_3.14.0          
##  [3] CoreGx_2.7.4                SummarizedExperiment_1.36.0
##  [5] Biobase_2.66.0              GenomicRanges_1.58.0       
##  [7] GenomeInfoDb_1.42.0         IRanges_2.40.0             
##  [9] S4Vectors_0.44.0            MatrixGenerics_1.18.0      
## [11] matrixStats_1.5.0           BiocGenerics_0.52.0        
## [13] UpSetR_1.4.1                scales_1.4.0               
## [15] ggplot2_4.0.3               data.table_1.18.4          
## 
## loaded via a namespace (and not attached):
##  [1] bitops_1.0-9            gridExtra_2.3.1         rlang_1.3.0            
##  [4] magrittr_2.0.5          shinydashboard_0.7.3    otel_0.2.0             
##  [7] compiler_4.4.3          reshape2_1.4.5          vctrs_0.7.3            
## [10] relations_0.6-17        stringr_1.6.0           pkgconfig_2.0.3        
## [13] crayon_1.5.3            fastmap_1.2.0           backports_1.5.1        
## [16] XVector_0.46.0          labeling_0.4.3          caTools_1.18.3         
## [19] promises_1.5.0          rmarkdown_2.31          UCSC.utils_1.2.0       
## [22] coop_0.6-3              xfun_0.59               zlibbioc_1.52.0        
## [25] cachem_1.1.0            jsonlite_2.0.0          SnowballC_0.7.1        
## [28] later_1.4.8             DelayedArray_0.32.0     BiocParallel_1.40.0    
## [31] parallel_4.4.3          sets_1.0-25             cluster_2.1.8.2        
## [34] R6_2.6.1                stringi_1.8.7           bslib_0.11.0           
## [37] RColorBrewer_1.1-3      limma_3.62.2            boot_1.3-32            
## [40] jquerylib_0.1.4         Rcpp_1.1.2              bookdown_0.47          
## [43] knitr_1.51              downloader_0.4.1        httpuv_1.6.17          
## [46] Matrix_1.7-5            igraph_2.3.3            tidyselect_1.2.1       
## [49] abind_1.4-8             yaml_2.3.12             gplots_3.3.0           
## [52] codetools_0.2-20        lattice_0.22-9          tibble_3.3.1           
## [55] plyr_1.8.9              shiny_1.14.0            withr_3.0.3            
## [58] BumpyMatrix_1.10.0      S7_0.2.2                evaluate_1.0.5         
## [61] bench_1.1.4             pillar_1.11.1           lsa_0.73.4             
## [64] KernSmooth_2.23-26      checkmate_2.3.4         DT_0.34.0              
## [67] shinyjs_2.1.1           piano_2.22.0            generics_0.1.4         
## [70] gtools_3.9.5            xtable_1.8-8            marray_1.84.0          
## [73] glue_1.8.1              slam_0.1-55             tools_4.4.3            
## [76] fgsea_1.32.4            visNetwork_2.1.4        fastmatch_1.1-8        
## [79] cowplot_1.2.0           grid_4.4.3              GenomeInfoDbData_1.2.13
## [82] cli_3.6.6               viridisLite_0.4.3       S4Arrays_1.6.0         
## [85] dplyr_1.2.1             gtable_0.3.6            sass_0.4.10            
## [88] digest_0.6.39           SparseArray_1.6.0       htmlwidgets_1.6.4      
## [91] farver_2.1.2            htmltools_0.5.9         lifecycle_1.0.5        
## [94] httr_1.4.8              statmod_1.5.2           mime_0.13