Level 0

The “mc0” and “sc0” tables are identical, other than containing m0id rather than s0id, respectively. See the Level 0 Pre-processing section for more information.

Information about the different well types is available in Level 0 Pre-processing section.

SC Data-containing Tables

- Level 1

- Level 2

- Aggregated IDs

- Representative Samples

See the Data Interpretation>Representative Samples section for more details.

MC Data-containing Tables

- Level 1

- Level 2

- Level 3

- Aggregated IDs

- Level 4

- Level 4 Parameters

- Level 5

See the Data Interpretation>Hit Calls section for more details.

- Level 5 Parameters

- Representative Samples

See the Data Interpretation>Representative Samples section for more details.

- Level 6

See the Data Interpretation>Flags section for more details.

- Level 7

See the Data Interpretation>Adminstered Equivalent Doses section for more details.

Assay and Auxiliary Annotation Tables

The fields pertinent to the tcpl package are listed in the tables below. More specifics on assay and auxiliary annotations will be provided in later sections.

** indicates tables may have limited tcpl functionality, but data is still retrievable via tcplQuery.

- Assay Source

- Assay

- Assay Component

- Assay Component Endpoint

See the Data Interpretation>Cytotoxicity Burst Distribution section for more details on “burst_assay”.

- Assay Component Map

- Assay Descriptions

See the Data Interpretation>Assay Description Documents section for more details.

- Chemical

Chemical ID is the DSSTox GSID within the ToxCast data, but can be any integer and will be auto-generated (if not explicitly defined) for newly registered chemicals.

- Chemical Lists

- Sample

The stock concentration fields in the “sample” table allow the user to track the original concentration when the neat sample is solubilized in vehicle before any serial dilutions for testing purposes. The US EPA’s ChemTrack application and database supports chemical procurement and sample management for ToxCast in vitro screening efforts.

Assay Nomenclature

The definition of an “assay” is, for the purposes of this package, broken into:

• assay_source: the vendor/origination of the data
  
• assay: the procedure to generate the component data
• assay_component: the raw data readout(s)
  
• assay_component_endpoint: the normalized component data

Assay source, assay, assay component, and assay endpoint are registered via tcpl scripting into a collection of tables. These assay element tables broadly describe who conducted the assay, what platform was used, what was being measured (raw readout), and how the measurement was interpreted (normalized component data). A hierarchical structure of the assay elements is as follows: assay source > assay > assay component > assay component endpoint.

As one moves down the hierarchy, each additional level has a ‘one-to-many’ relationship with the previous level. For example, an assay component can have multiple assay endpoints, but an assay endpoint can derive only from a single assay component.

Minimum Required Fields

Throughout the tcpl R package, the levels of assay hierarchy are defined and referenced by their auto-incremented primary keys in the tcpl database: asid (assay source ID), aid (assay ID), acid (assay component ID), and aeid (assay endpoint ID). These abbreviations mirror the abbreviations for the identifiers (ids) with “nm” in place of “id” in the abbreviations, e.g. assay_component_name is abbreviated acnm.

All processing occurs by assay component or assay endpoint, depending on the processing type (single-concentration or multiple-concentration) and level. No data is stored at the assay or assay source level. The “assay” and “assay_source” tables store annotations to help in the processing and down-stream understanding of the data. Additional details for registering each assay element and updating annotations are provided below. In addition to each assay element’s id, the minimal registration fields in order to ‘pipeline’ are:

• assay_source_name (asnm)
• assay_name (anm)
• assay_footprint
• assay_component_name (aenm)
• assay_component_endpoint_name (aenm)
• normalized_data_type

Assay Source

Assay source refers to the vendor or origination of the data. To register an assay source, an unused asid must be selected to prevent overwriting of existing data. When adding a new assay source, this should be an abbreviation, as subsequent levels will build on this assay source name.

tcplLoadAsid()
tcplRegister(what = "asid", flds = list(asid = 1, asnm = "Tox21"))

The tcplRegister function takes the abbreviation for assay_source_name, but the function will also take the unabbreviated form. The same is true of the tcplLoadA- functions, which load the information for the assay annotations stored in the database.

Assay

Assay refers to the procedure, conducted by some vendor, to generate the component data. To register an assay, an asid must be provided to map the assay to the correct assay source. One source may have many assays. To ensure consistency of the naming convention, first check how other registered assays within the assay source were conducted and named. The assay names follow an abbreviated and flexible naming convention of Source_Assay. Notable assay design features to describe the assay include:

• Technology (i.e., detection technology),
• Format (e.g., organism, tissue, cell short name, or cell free component source name),
• Target (i.e., intended target, intended target family, gene), or
• Objective aspects (e.g., timepoint or assay footprint).

The most distinguishing features will be selected to create a succinct assay name. Variation depends on the assay itself as well as other assays provided by the vendor. If multiple features are needed to describe an assay, order will be based on relative importance in describing the assay and the assay’s relation to other assays provided by the vendor to limit confusion. “Source_Technology_Format_Target” is a commonly used naming order. However, if one target is screened on different assay platforms by the vendor, “Source_Target_Technology” is a more appropriate naming convention. This is the case for the Tox21 assays. Additional features may be relevant, including agonist or antagonist mode, or “Follow-up” if the assay is a secondary specificity assay. Conversely, some assays utilize a cell-based format to screen a functional profile of targets. These assays follow a naming convention, Source_Format, where specific target information is defined at the component and endpoint level. Bioseek and Attagene are sources that provide cell-based assays. Considering the diversity of the data sources and high throughput assays in ToxCast, a flexible naming approach is best used in conjunction with subject matter expert discretion.

tcplLoadAid(what = "asid", val = 1)
tcplRegister(what = "aid", flds = list(asid = 1, anm = "TOX21_ERa_BLA_Agonist", assay_footprint = "1536 well"))

When registering an assay (aid), the user must give an asid to map the assay to the correct assay source. Registering an assay, in addition to an assay_name (anm) and asid, requires assay_footprint. The assay_footprint field is used in the assay plate visualization functions (discussed later) to define the appropriate plate size. The assay_footprint field can take most string values, but only the numeric value will be extracted, e.g. the text string “hello 384” would indicate to draw a 384-well microtitier plate. Values containing multiple numeric values in assay_footprint may cause errors in plotting plate diagrams.

Assay Component

Assay component, or “component” for short, describes the raw data readouts. Like the previous level, one assay may have many components. To register an assay component and create an acid, an aid must be provided to map the component to the correct assay. The assay component name will build on its respective assay name, to describe the specific feature being measured in each component. If there is only one component, the component name can be the same as the assay name. If there are multiple components measured in an assay, understanding the differences, and how one component may relate to another within an assay, are important naming considerations to prevent confusion. Assay component names will usually follow the naming convention of Source_Assay_Component, where “Component” is a brief description of what is being measured.

tcplLoadAcid(what = "asid", val = 1, add.fld = c("aid", "anm"))
tcplRegister(what = "acid", flds = list(aid = 1, acnm = "TOX21_ERa_BLA_Agonist_ratio"))

The final piece of assay information needed is the assay component source name (abbreviated acsn), stored in the “assay_component_map” table. The assay component source name is intended to simplify level 0 pre-processing by defining unique character strings (concatenating information if necessary) from the source files that identify the specific assay components. An assay component can have multiple acsn values, but an acsn must be unique to one assay component. Assay components can have multiple acsn values to minimize the amount of data manipulation required (and therefore potential errors) during the level 0 pre-processing if assay source files change or are inconsistent. The unique character strings (acsn) get mapped to acid.

tcplRegister(what = "acsn", flds = list(acid = 1, acsn = "TCPL-MC-Demo"))

Assay Component Endpoint

Assay component endpoint, or “endpoint” for short, represents the normalized component data. To register an endpoint and create an aeid, an acid must be provided to map the endpoint to the correct component. In past tcpl versions, each component could have up to two endpoints therefore endpoint names would express directionality (*_up/_down*). tcpl v3+ allows bidirectional fitting to capture both the gain and loss of signal. Therefore with tcpl v3+ , the endpoint name will usually be the same as the component name.

tcplLoadAeid(fld = "asid", val = 1, add.fld = c("aid", "anm", "acid", "acnm"))
tcplRegister(what = "aeid", flds = list(acid = 1, aenm = "TOX21_ERa_BLA_Agonist_ratio", normalized_data_type = "percent_activity", export_ready = 1, burst_assay = 0))

Registering an assay endpoint also requires the normalized_data_type field. The normalized_data_type is used when plotting and currently, the package supports the following values: percent_activity, log2_fold_induction, log10_fold_induction, and fold_induction. Any other values will be treated as “percent_activity.”

Other required fields to register an assay endpoint do not have to be explicitly defined and will default to 0 if not provided. These fields represent Boolean values (1 or 0, 1 being TRUE ). The export_ready field indicates (1) the data is done and ready for export or (0) still in progress. The burst_assay field is specific to multiple-concentration processing and indicates (1) the assay endpoint is included in the burst distribution calculation or (0) not.

Naming Revision

There are circumstances where assay, assay component, and assay endpoint names change. The aid, acid, and aeid are considered more stable in the database, and these auto-incremented keys should not change. To revise naming for assay elements, the correct id must be specified in the tcplUpdate statement to prevent overwriting data.

tcplUpdate(what = "acid", flds = list(aid = 1, acnm = "TOX21_ERa_BLA_Agonist_ratio"))

Reasons for name changes could include:

• Feedback from subject matter experts or assay data generators;
• Clarifications on cell line and cell line drift;
• Addition of new assay data that makes the old naming convention insufficient, such as antagonist assays run with different concentrations of an agonist; or
• Laboratory or Center reorganizations.

Thus, users should be advised that while assay naming is used to infer information about the biology of the assay, assay naming will change over time to reflect progress in building ToxCast as a data resource.

Chemical Lists

Optionally, the user can subdivide the chemical IDs based on presence in different chemical lists using tcplLoadChemList. These chemical lists are curated by the US EPA in the Distributed Structure-Searchable Toxicity (DSSTox) database. Chemicals can belong to more than one chemical list, and will be listed as separate entries when loading chemical list information.

tcplLoadChemList(field = "chid", val = 1:2)

Required Fields

The N/A column indicates whether the field can be NULL in the pre-processed data when writing Level 0. In past versions of tcpl, there were some exceptions where concentrations could be NULL. All conc values must be numeric for processing. For blank and neutral control wells, NULL concs will automatically be set to zero if not provided. If other concs are NULL, user will get a warning and be unable to write Level 0. Additionally, if the raw value is NULL, well quality must be 0.

The well type (wllt) field is used to differentiate wells in numerous applications, including normalization and definition of the assay noise level. Package users are encouraged to utilize additional well types (for example, well types “x”, “y”, “z”) or suggest new methods to better accommodate their data.

Well Types

Before writing and processing any data to the database, the user must register the assay and chemical information, as described above. tcpl includes three functions for adding new data:

• tcplRegister: to register a new assay element or chemical
• tcplUpdate: to change or add additional information for existing assay or chemical ids
• tcplWriteLvl0: to load formatted source Level 0 data

The final step in level 0 pre-processing is to load the data into the database. The tcplWriteLvl0 function loads data into the database by checking each required field for expected input, such as correct class or registered sample IDs and concentration values for test wells. Assay component source name, if provided, will map to the appropriate acid. The type argument is used throughout the package to distinguish the screening paradigm (SC or MC) and therefore processing required.

# Write/load the 'mcdat' into the database.
tcplWriteLvl0(dat = mcdat, type = "mc")

As a check to confirm data has been written, the tcplLoadData function is used to load data from the database into user’s R session. Furthermore, the tcplPrepOtpt function can be used in conjunction with tcplLoadData to prepare the data in a readable format with additional chemical and assay annotation information. See Data Retrieval sections for further details.

# Load the level 0 data from the database to R.
tcplLoadData(lvl = 0, fld = "acid", val = 1, type = "mc")
tcplPrepOtpt(tcplLoadData(lvl = 0, fld = "acid", val = 1, type = "mc"))

In the loaded Level 0 data, acsn is replaced with correct acid and the m0id field is added. These “m#” fields in MC data are primary keys for each level of data and can link the various levels of data. All keys are auto-generated and will change anytime data are reprocessed. Note, the primary keys only change for the levels affected, e.g. if the user reprocesses level 1, the level 0 keys will remain the same.

Overview

All processing in the tcpl package occurs at the assay component or assay endpoint level. There is no capability in either SC or MC processing to combine data from multiple assay components or assay endpoints. Any combining of data must occur before or after processing. For example, a ratio of two raw values could be processed if the user calculated the ratio during the custom pre-processing and uploaded values as a single “component”.

Once the Level 0 data are loaded, data processing occurs via the tcplRun function for both SC and MC multiple-concentration screening. tcplRun can either take a single ID (acid or aeid, depending on the processing type and level) or an asid. If given an asid, the tcplRun function will attempt to process all corresponding components/endpoints. When processing by acid or aeid, the user must supply correct ID for each level.

The processing is sequential, and every level of processing requires successful processing at the antecedent level. Any processing changes will trigger a “delete cascade,” removing any subsequent data affected by the processing change to ensure complete data fidelity. For example, processing level 3 data will first require data from levels 4 through 6 to be deleted for the corresponding IDs. Changing method assignments will also trigger a delete cascade for any corresponding data.

For tcplRun, the user must supply a starting level (slvl) and ending level (elvl). There are four phases of processing, as reflected by messages printed in console: (1) data for the given IDs are loaded, (2) the data are processed, (3) data for the same ID in subsequent levels are deleted, and (4) the processed data is written to the database. The ‘outfile’ parameter can give the user the option of print this output text to a file. If an id fails while processing multiple levels, the function will not attempt to process the failed ids in subsequent levels. When finished processing, a list indicating the processing success of each ID is returned. For each level processed, the list will contain two elements: (1) “l#” a named Boolean vector where TRUE indicates successful processing, and (2) “l#_failed” containing the names of any ids that failed processing, where “#” is the processing level.

Processing of multiple assay components or endpoints can be executed simultaneously. This is done with the internal utilization of the mclapply function from the parallel package. Parallel processing is done by id. Depending on the system environment and memory constraints, the user may wish to use more or less processing power. For processing on a Windows operating system, the default is mc.cores = 1, unless otherwise specified. For processing on a Unix-based operating system, the default is mc.cores = NULL i.e. to utilize all cores except for one, which is necessary for ‘overhead’ processes. The user can specify more or less processing power by setting the “mc.cores” parameter to the desired level. Note, this specification should meet the following criteria $1 \space \leq \space \mathit{mc.cores} \space \leq \space \mathit{detectCores()}-1$.

Level IDs

In this table, the Input ID column indicates the ID used for each processing step where Method ID indicates the ID used for assigning methods for data processing, when necessary. SC1 requires an acid input, but the methods are assigned by aeid. The same is true for MC3 processing. SC1 and MC3 are the normalization steps and convert acid to aeid. Only MC2 has methods assigned by acid.

The processing requirements vary by screening paradigm and level, which later sections cover in detail. However, in general, specific method assignments will be required to accommodate different experimental designs or data processing approaches.

Methods

To promote reproducibility, all method assignments are saved in the database and utilize methods described in the available list of methods for each processing level. In general, methods data are stored in the “_methods” and “_id” tables for each corresponding level. For example, the “sc1” table is accompanied by the “sc1_methods” table which stores the available list of methods for SC1, and the “sc1_aeid” table which stores the method assignments and execution order.

There are three functions to easily modify and load method assignments:

• tcplMthdAssign – assigns methods to specified id(s)
• tcplMthdClear – clears method assignments to specified id(s)
• tcplMthdLoad – queries database and returns the method assignments for specified id(s)
  
• tcplMthdList – queries database and returns available methods at specified level(s)
  

The following code blocks provides examples of the method-related functions:

## Methods Assignment ##
# For illustrative purposes, assign level 2 MC methods to ACIDs 97, 98, and 99.
# First check for available methods.
mthds <- tcplMthdList(lvl = 2, type = "mc")
mthds[1:2]
# Assign some methods to ACID 97, 98, and 99.
tcplMthdAssign(lvl = 2, 
               id = 97:99, 
               mthd_id = c(3, 4, 2), 
               ordr = 1:3, 
               type = "mc")
# Check the assigned methods for ACID 97, 98, and 99 in the database.
tcplMthdLoad(lvl = 2, id = 97:99, type = "mc")

# Methods can be cleared one at a time for the given id(s)
tcplMthdClear(lvl = 2, id = 99, mthd_id = 2, type = "mc")
# Check the assigned methods for the single id updated, namely ACID 99.
tcplMthdLoad(lvl = 2, id = 99, type = "mc")

# Or all methods can be cleared for the given id(s)
tcplMthdClear(lvl = 2, id = 97:98, type = "mc")
# Check the assigned methods for the all updated ids, namely ACID 97 and 98.
tcplMthdLoad(lvl = 2, id = 97:98, type = "mc")

Later sections of this vignette provide level-specific method assignment examples and more details on the methods themselves. Most examples will reflect commonly-used methods assigned, but one should consider the data at hand and all methods available for the level prior to assigning.

Data Normalization

Data normalization occurs in both SC and MC processing paradigms at levels 1 and 3, respectively. While the two paradigms use different methods, the normalization approach is the same for both. Data normalization does not have to occur within the package as pre-normalized data can be loaded into the database at Level 0. However, data must be zero-centered. Thus, the data must either be zero-centered in Level 0 pre-processing or the user must pick a methodology from the associated level 1 and 3 methods to zero-center the data before model fitting occurs.

Fold-change and a percent of control are typical approaches to normalization. Given data must be zero-centered, fold-change data in general is log-transformed. Log-scale transformations for fold-change data is typically base 2 (log₂), but other bases may be more appropriate in some circumstances.

Normalizing to a percent of control requires three normalization methods:

1. one to define the baseline value,
2. one to define the control value, and
3. one to calculate percent of control (“resp.pc”).

Normalizing to fold-change also requires three methods: 1. one to define the baseline value, 2. one to calculate the fold-change, and 3. one to log-transform the fold-change values (“resp.fc”).

(cval) is the corrected response value for a test well defined in MC2. Methods defining a baseline value (bval) have the “bval” prefix, methods defining the positive control value (pval) have the “pval” prefix. Pval may be set as 0 if no positive control wells are provided and measuring decreases in signal. Finally, methods that calculate the final response value have the “resp” prefix. For example, “resp.log2” does a log-transformation of the response value using a base value of 2. The formulae for calculating the percent of control and fold-change response are listed in equations 1 and 2, respectively. Note that the fold-change calculation divides by the baseline value and thus must have some non-zero values associated with the baseline to successfully calculate fold-change.

$$ resp.pc = \frac{cval - bval}{pval - bval}*100 $$

$$ resp.fc = \frac{cval}{bval} $$

Order matters when assigning normalization methods. The bval, and pval if normalizing as a percent of control, need to be calculated prior to calculating the response value. Examples of normalization schemes are presented below:

If the data does not require any normalization, the “none” method will be assigned. The “none” method simply copies the input data to the response field. Without assigning “none”, the response field will not get generated and processing will fail.

With tcpl v2 , responses were only fit in the positive analysis direction. Therefore, a signal in the negative direction needed to be “flipped” to the positive direction during normalization. Multiple endpoints stemming from one component were created to enable multiple normalization approaches when the assay measured gain and loss of signal. Negative direction data was inverted by multiplying the final response values by −1 via the “resp.multneg1” methods. For tcpl v3 onward, the tcplFit2 package is utilized which allows for bidirectional fitting, meaning the “resp.multneg1” method is now only required in special cases.

In addition to the required normalization methods, the user can apply additional methods to transform the normalized values. For example, “resp.blineshift.50.spid” corrects for baseline deviations by spid. A complete list of available methods, by processing type and level, can be accessed with tcplMthdList. More information is also available in package documentation, ??tcpl::Methods.

SC Data Processing

The goal of single-concentration processing is to identify potentially active compounds from a broad screen at a single concentration. SC processing consists of 2 levels:

> Level 1

Level 1 processing converts the assay component to assay endpoint(s), defines the normalized-response value (resp), and optionally, derives the baseline value (bval) and positive control value (pval). The purpose of level 1 is to normalize the raw values to either the percentage of a control or fold-change from baseline.

- Methods Assignment

The first step in beginning the processing is to identify which assay endpoints stem from the assay component(s) being processed. With the corresponding endpoints identified, the appropriate methods can be assigned.

# Load the 'aeid' values for acid 2 ##
tcplLoadAeid(fld = "acid", val = 2)

# Assign the level 1 methods to aeid 1 and 2 ##
tcplMthdAssign(lvl = 1,  # processing level
               id = 1:2, # assay endpoint ID's to assign methods
               mthd_id = c(1, 11, 13), # method(s) to be assigned
               ordr = 1:3, # order the method(s) should be applied
               type = "sc") # the data/processing type

Above, methods 1, 11, and 13 were assigned for both endpoints. The method assignments instruct the processing to: (1) calculate bval for each assay plate ID by taking the median of all data where the well type equals “n”; (2) calculate a fold-change with respect to the bval (i.e. $\frac{resp}{\mathit{bval}}$); (3) log-transform the fold-change values with base 2. For a complete list of normalization methods see tcplMthdList(lvl = 1, type = “sc”) or ?SC1_Methods .

If a user needs to add a method to the end of a normalization sequence, as shown above, then the user can use a second method assignment statement. For example, if AEID 2 should indicate a change of directionality to be more biologically interpretable, all responses can be multiplied by −1. Reminder, the order of methods assignment matters, particularly in the normalization step.

# Assign a fourth step to the normalization processing - for AEID 2 only.
tcplMthdAssign(lvl = 1, # processing level
               id = 2, # assay endpoint ID's to assign methods
               mthd_id = 16, # method(s) to be assigned
               ordr = 4, # order the method(s) should be applied
               type = "sc") # the data/processing type

With the normalization methods defined, the data are ready for SC1 processing. Before normalization occurs within tcplRun, all wells with well quality (wllq) 0 are removed.

> Level 2

Level 2 processing defines the baseline median absolute deviation (bmad), collapses any replicates by sample ID (spid), and determines the activity.

- Methods Assignment

Before the data are collapsed by spid, the bmad is calculated as the median absolute deviation of all treatment wells (wllt = t - default option) or neutral control wells (wllt = n). The calculation defines bmad for the entire endpoint. If additional data is added, the bmad for the associated endpoints will be recalculated. Note, this bmad equation is different from MC screening.

bmad_sc = 1.4826 * median(|y_i − ỹ|)

Where y_i is the i^th observation of all wells within a well type in the assay component and ỹ is the median across all y_i’s. The constant value, 1.4826, is the default adjustment value used in the underlying R function to ensure bmad is a consistent estimator of the standard deviation (σ) assuming the sample size (N) of the observations is large and they are normally distributed (i.e. Gaussian), see mad() in R and unbiased mad for further details.

To collapse the data by spid, the median response of replicates at each concentration index is calculated. SC data screening involves testing at 1-3 concentrations, which is insufficient for fitting. Data are then further collapsed by taking the maximum of those median values (max_med). Once collapsed, such that each endpoint-sample has one value, the activity is determined. For an active hit call (hitc), the sample’s max_med must be greater than a specified efficacy cutoff (coff).

In the below code, methods are assigned such that the cutoff value is log₂(1.2). Thus, if the maximum median value (max_med) is greater than or equal to the efficacy cutoff (coff = log₂(1.2)), then the sample ID is considered active and the hit call (hitc) is set to 1. The coff is defined as the maximum of all values given by the assigned level 2 methods. Failing to assign a level 2 method will result in every sample being called active. For a complete list of level 2 methods, see tcplMthdList(lvl = 2, type = “sc”) or ?SC2_Methods. With the methods assigned and the cutoff set, the data are ready for SC2 processing.

# Assign a cutoff value of log2(1.2) ##
tcplMthdAssign(lvl = 2, # processing level
               id = 1,  # assay endpoint ID's to assign methods
               mthd_id = 3, # method(s) to be assigned 
               type = "sc") # the data/processing type

## SC2 processing for acid 1 ##
sc2_res <- tcplRun(id = 1, slvl = 2, elvl = 2, type = "sc")

Sometimes it is necessary to prevent hitcalling for responses in the biologically irrelevant direction. SC2 contains two methods for overwriting the max_med value. If applied, negative hitcalls will be given for any max_med greater or less than cutoff in the biologically unintended direction, by comparing the max_med to either the positive or negative cutoff instead examining absolute values.

MC Data Processing

The goal of multiple-concentration processing is to estimate the activity, potency, efficacy, and other parameters for sample-assay pairs.

> Level 1

Level 1 processing defines the replicate and concentration index fields to facilitate downstream processing. Because of cost, availability, physicochemical, and technical constraints, screening efforts may utilize different experimental designs. The resulting data may contain an inconsistent number of concentration groups, concentration values, and technical replicates. To enable quick and uniform processing, level 1 processing explicitly defines concentration and replicate indices as 1…N for increasing concentrations and technical replicates, where 1 represents the lowest concentration or first technical replicate.

To assign replicate and concentration indices, we assume one of two experimental designs. The first design assumes samples are plated in multiple concentrations on each assay plate, such that the concentration series all fall on a single assay plate. The second design assumes samples are plated in a single concentration on each assay plate, such that the concentration series falls across many assay plates.

For both experimental designs, data are ordered by source file (srcf), assay plate ID (apid), column index (coli), row index (rowi), sample ID (spid), and concentration (conc). Concentration is rounded to three significant figures to correct for potential rounding errors. After ordering the data, a temporary replicate ID is created for each concentration series. For test compounds in experimental designs with the concentration series on a single plate and all control compounds, the temporary replicate ID consists of the sample ID, well type (wllt), source file, assay plate ID, and concentration. The temporary replicate ID for test compounds in experimental designs with concentration series that span multiple assay plates is defined similarly, but does not include the assay plate ID.

Once the data are ordered, and the temporary replicate ID is defined, the data are scanned from top to bottom and the replicate index (repi) incremented every time a replicate ID is duplicated. Then, for each replicate, the concentration index (cndx) is defined by ranking the unique concentrations, with the lowest concentration starting at 1. No methods need to be applied and the following demonstrates how to carry out the MC1 processing and look at the resulting data:

The package also contains a function, tcplPlotPlate , for visualizing the data at the assay plate level. This function can be used to visualize the data at levels 1 to 3.

tcplPlotPlate(dat = m1dat, apid = "4009721")

Design of Assay Plate 4009721

In the generated figure, the row and column indices are printed along the respective edges of the plate, with the raw observed values in each well represented by color. While the plate does not give sample ID information, the letter/number codes in the wells indicate the well type and concentration index, respectively. Wells with poor well quality (wllq==0 in Level 0) will display with an “X.” The title of the plate display lists the assay component/assay endpoint and the assay plate ID (apid).

> Level 2

Level 2 processing removes data where the well quality (wllq) equals 0 and defines the corrected value (cval) field. MC2 also allows for additional transformation of the raw values at the assay component level. Examples of transformations include basic arithmetic manipulations to complex spatial noise reduction algorithms, such as aggregation across biological replicates.

- Methods Assignment

Every assay component needs at least one transformation method assigned to complete level 2 processing, even if no transformations are necessary. In the following example, the “none” method will be assigned so MC2 processing can be completed.

## Methods Assignment ##
# Assign the level 2 transformation method 'none' to ACID 1.
tcplMthdAssign(lvl = 2, # processing level
               id = 1, # assay component ID's to assign methods
               mthd_id = 1, # method(s) to be assigned
               ordr = 1, # order of the method(s) should be assigned
               type = "mc") # the data/processing type

## MC2 processing for acid 1 ##
mc2_res <- tcplRun(id = 1, slvl = 2, elvl = 2, type = "mc")

For the complete list of level 2 transformation methods currently available, see tcplMthdList(lvl = 2, type = “mc”) or ?MC2_Methods for more details. The coding methodology used to implement the methods is beyond the scope of this vignette, but, in brief, the method names in the database correspond to a function name in the list of functions returned by mc2_mthds() (the mc2_mthds function is not exported, and not intended for use by the user). Each of the functions in the list given by mc2_mthds only return expression objects that the processing function called by tcplRun executes in the local function environment to avoid making additional copies of the data in memory. We encourage suggestions for new methods.

> Level 3

Level 3 processing converts the assay component to assay endpoint(s) and defines the normalized-response value field (resp); and optionally, the baseline value (bval) and positive control value (pval) fields. MC3 processing normalizes the corrected values to either the percentage of a control or to fold-change from baseline. The normalization process is discussed in greater detail in the Data Normalization section. A primary distinction between MC2 and MC3 processing is ID for each level.

- Methods Assignment

The user first needs to check which assay endpoints stem from the the assay component queued for processing. With the corresponding aeids identified, the normalization methods can be assigned. In the following example, methods 17, 9, and 7 were assigned for both endpoints. These methods involve: (1) calculating bval for each assay plate ID by taking the median of all data where the well type equals “n” or the well type equals “t” and the concentration index is 1 or 2; (2) calculating a fold-change over bval; (3) log-transforming the fold-change values with base 2. For a complete list of normalization methods see tcplMthdList(lvl = 3, type = “mc”) or ?MC3_Methods . With normalization methods defined, the data are ready for MC3 processing.

# Look at the assay endpoints for acid 1 ##
tcplLoadAeid(fld = "acid", val = 1)

## Methods Assignment ##
# Assign the baseline calculation and normalization methods to aeids 1 and 2.
tcplMthdAssign(lvl = 3, # processing level
               id = 1:2, # assay endpoint ID to assign methods
               mthd_id = c(17, 9, 7), # method(s) to be assigned
               ordr = 1:3, # order the method(s) should be applied
               type = "mc") # the data/processing type

## MC3 processing for acid 1 ##
mc3_res <- tcplRun(id = 1, slvl = 3, elvl = 3, type = "mc")

Notice that MC3 processing takes an acid, not an aeid, as the input ID. As mentioned in previous sections, the user must will assign MC3 normalization methods by aeid then process by acid. The MC3 processing will attempt to process all endpoints for a given component. If one endpoint fails for any reason (e.g., does not have appropriate methods assigned), the processing for the entire component fails.

> Level 4

Level 4 processing models the activity of each concentration-response series. Each series is bidirectionally fit using methods available in the tcplfit2 R package (Sheffield et al., 2021). Bidirectional fitting means curve inversion is not necessary as in past versions (i.e. observed ‘negative/decreasing’ responses multiplied by −1).

– Pre-Modeling Processes

Level 4 processing establishes a noise-band for the endpoint using the baseline median absolute deviation (bmad). Here, the bmad is calculated by the baseline response values, either untreated control wells (e.g. ‘wllt = n’; neutral solvent wells like DMSO) or test samples from the two lowest concentrations (i.e. ‘wllt = t’ & concentration index is 1 or 2). The bmad calculation is done across the entire endpoint. If additional data is added, the bmad values for all associated assay endpoints is recalculated.

bmad_mc = 1.4826 * median(|y_i − ỹ|) Where y_i is the i^th baseline observation as defined by the assigned method and ỹ is the median of all the baseline observations. The constant value, 1.4826, is the default adjustment value used in the underlying R function to ensure bmad is a consistent estimator of the standard deviation (σ) assuming the sample size (N) of the baseline observations is large and they are normally distributed (i.e. Gaussian), see mad() in R and unbiased mad for more details.

One standard deviation of baseline response is also calculated, since it is necessary for the calculation of the benchmark response (BMR) in both the curve-fitting and hit-calling functions from tcplfit2 , tcplfit_core and tcplhit2_core , respectively. One standard deviation of the baseline response will be estimated using the two lowest concentration groups of test samples or neutral control wells across all chemicals, depending on bmad method selected. BMR, like bmad, may change if additional data is added. Thus, the BMD estimates may also change. See the Data Interpretation>Benchmark Dose section for more details.

Before the model parameters are estimated, a set of summary values are calculated for each concentration-response series:

• Minimum and maximum observed responses (resp_min & resp_max, respectively).
• Minimum and maximum concentrations (conc_min & conc_max, respectively).
• The total number of concentration groups (nconc).
• Total number of observed responses (i.e. data points in the concentration series) (npts).
• Number of replicates in concentration groups (nrep).
• The maximum mean and median responses along with the concentration at which they occur (max_mean, max_med, max_mean_conc, & max_med_conc, respectively).
• The minimum mean and median responses along with the concentration at which they occur (min_mean, min_med, min_mean_conc, & min_med_conc, respectively).
• The maximum median difference response – the greater distance to 0 of max_med and min_med – along with the concentration at which it occurs (max_med_diff & max_med_diff_conc, respectively)
• The number of median responses greater than 3 * bmad or less than −3 * bmad (nmed_gtbl_pos & nmed_gtbl_neg, respectively).

The “mean” and “median” response values are defined as the mean or median of all observed response values at each concentration. In other words, the maximum median is the maximum of all median response values across all concentrations in the series. Similarly, the minimum median is the true minimum (or max in the negative direction) of all median response values across the concentrations in the series.

The following code demonstrates how to load the MC3 data for a single aeid.

## Evaluate the MC3 Data ##
# Load the MC3 data from the database.
mc3 <- tcplLoadData(lvl = 3,
                    type = 'mc',
                    fld = 'aeid',
                    val = 80)
# Prepare the data into a readable format.
mc3 <- tcplPrepOtpt(mc3)

For demonstration purposes, the mc_vignette R data object is provided in the package since the vignette is not directly connected to such a database. The mc_vignette object contains a subset of data from levels 3 through 5 from invitrodb v4.2. The following code loads the example mc3 data object, then plots the concentration-response series for an example spid with the summary estimates indicated.

# Load the example data from the `tcpl` package.
data(mc_vignette, package = 'tcpl')
# Allocate the level 3 example data to `mc3`.
mc3_example <- mc_vignette[['mc3']]
# level 3 does not store logc anymore, create it for plotting purposes
mc3_example[, logc := log10(conc)]
# Obtain the MC4 example data.
mc4_example <- mc_vignette[["mc4"]]
# Obtain the minimum response observed and the 'logc' group - 'resp_min'.
level3_min <- mc3_example %>%
  dplyr::group_by(spid, chnm) %>% 
  dplyr::filter(resp == min(resp)) %>% 
  dplyr::filter(spid == "01504209")
# Obtain the maximum response observed and the 'logc' group - 'resp_max'.
level3_max <- mc3_example %>% 
  dplyr::group_by(spid, chnm) %>% 
  dplyr::filter(resp == max(resp)) %>% 
  dplyr::filter(spid == "01504209")
# Obtain the level 3 data and 'center' estimates for responses per 'logc' group.
level3_summary <- mc3_example %>%
  dplyr::filter(spid == "01504209") %>% 
  dplyr::select(., c(spid, chnm, logc, resp)) %>%
  dplyr::group_by(spid, chnm, logc) %>% 
  dplyr::summarise(mean_resp = mean(resp), med_resp = median(resp))

## Generate Individual Summary Plots ##
# Plot the mean responses for each log-concentration group.
A <- mc3_example %>%
  dplyr::filter(spid == "01504209") %>% 
  ggplot(data = ., aes(logc, resp)) +
  geom_point(pch = 1, size = 2) +
  geom_point(data = level3_summary,
             aes(x = logc, y = mean_resp,
                 col = 'mean responses'),
             alpha = 0.75,size = 2) +
  scale_color_manual(values = 'paleturquoise3',
                     aesthetics = 'col') +
  labs(lty = "", colour = "")+
  xlab(expression(paste(log[10],"(Concentration) ", mu, "M"))) +
  ylab(expression(paste(log[2], "(Fold Induction)"))) +
  ggtitle("Mean Responses") +
  theme_bw() +
  theme(legend.position = 'bottom')
# Plot the median responses for each log-concentration group.
B <- mc3_example %>%
  dplyr::filter(spid == "01504209") %>% 
  ggplot(data = .,aes(logc,resp)) +
  geom_point(pch = 1, size = 2) +
  geom_point(data = level3_summary,
             aes(x = logc, y = med_resp,
                 col = 'median response'),
             alpha = 0.75, size = 2) +
  scale_color_manual(values = 'hotpink',
                     aesthetics = 'col') +
  labs(lty = "", colour = "")+
  xlab(expression(paste(log[10], "(Concentration) ", mu, "M"))) +
  ylab(expression(paste(log[2], "(Fold Induction)"))) +
  ggtitle("Median Responses") +
  theme_bw() +
  theme(legend.position = 'bottom')
# Plot the maximum mean & median responses at the related log-concentration -
#   'max_mean' & 'max_mean_conc'.
C <- mc3_example %>%
  dplyr::filter(spid == "01504209") %>% 
  ggplot(data = .,aes(logc, resp)) +
  geom_point(pch = 1, size = 2) +
  geom_point(data = dplyr::filter(mc4, spid == "01504209"),
             aes(x = log10(max_mean_conc), y = max_mean,
                 col = 'maximum mean response'),
             alpha = 0.75, size = 2)+
  scale_color_manual(values = 'paleturquoise3',
                     aesthetics = 'col') +
  labs(lty = "", colour = "")+
  xlab(expression(paste(log[10], "(Concentration) ", mu, "M"))) +
  ylab(expression(paste(log[2], "(Fold Induction)"))) +
  ggtitle(label = "Maximum Mean Response") +
  theme_bw() +
  theme(legend.position = 'bottom')
# Plot the maximum mean & median responses at the related log-concentration -
#   'max_med' & 'max_med_conc'.
D <- example %>%
  dplyr::filter(spid == "01504209") %>% 
  ggplot(data = ., aes(logc, resp)) +
  geom_point(pch = 1, size = 2) +
  geom_point(data = dplyr::filter(mc4, spid == "01504209"),
             aes(x = log10(max_med_conc), y = max_med,
                 col = "maximum median response"),
             alpha = 0.75, size = 2)+
  scale_color_manual(values = 'hotpink',
                     aesthetics = 'col') +
  labs(lty = "", colour = "") +
  xlab(expression(paste(log[10], "(Concentration) ", mu, "M"))) +
  ylab(expression(paste(log[2], "(Fold Induction)"))) +
  ggtitle(label = "Maximum Median Response") +
  theme_bw() +
  theme(legend.position = 'bottom')
# Plot the minimum & maximum observed responses.
E <- mc3_example %>%
  dplyr::filter(spid == "01504209") %>% 
  ggplot(data = ., aes(logc, resp)) +
  geom_point(pch = 1, size = 2) +
  geom_point(data = level3_min,
             aes(x = logc, y = resp,
                 col = "minimum response"),
             alpha = 0.75, size = 2) +
  geom_point(data = level3_max,
             aes(x = logc, y = resp,
                 col = "maximum response"),
             alpha = 0.75, size = 2) +
  scale_color_manual(values = c('red', 'blue'),
                     aesthetics = 'col') +
  labs(lty = "", colour = "") +
  xlab(expression(paste(log[10], "(Concentration) ", mu,"M"))) +
  ylab(expression(paste(log[2], "(Fold Induction)"))) +
  ggtitle(label = "Minimum & Maximum\nResponses") +
  theme_bw() +
  theme(legend.position = 'bottom')
# Plot the minimum & maximum experimental log-concentration groups -
#   'logc_min' & 'logc_max'.
G <- mc3_example %>%
  dplyr::filter(spid == "01504209") %>% 
  ggplot(data = ., aes(logc, resp)) +
  geom_point(pch = 1, size = 2) +
  geom_vline(data = dplyr::filter(mc4, spid == "01504209"),
             aes(xintercept = log10(conc_min),
                 col = 'minimum concentration'),
             lty = "dashed") +
  geom_vline(data = dplyr::filter(mc4, spid == "01504209"),
             aes(xintercept = log10(conc_max),
                 col = 'maximum concentration'),
             lty = "dashed") +
  scale_color_manual(values = c('red', 'blue'),
                     aesthetics = 'col') +
  labs(lty = "", colour = "") +
  xlab(expression(paste(log[10], "(Concentration) ", mu, "M"))) +
  ylab(expression(paste(log[2], "(Fold Induction)"))) +
  ggtitle(label = "Minimum & Maximum\nConcentrations") +
  theme_bw() +
  theme(legend.position = 'bottom')
## Compile Summary Plots in One Figure ##
gridExtra::grid.arrange(
  A,B,C,D,E,G,
  nrow = 3, ncol = 2,
  top = mc3[which(mc4[,spid] == "01504209"), aenm]
)

These plots illustrate summary estimates calculated as part of the level 4 processing, which occurs prior to dose-response modeling. Each plot depicts the observed concentration-response data as white circles, where the x-axis is base 10 log-transformed concentration values. In the upper plots, the mean response values for each concentration group is depicted as turquoise circles (left) where as median response values for each concentration group are hot-pink circles (right). The middle plots depict mean and median responses, but only shows the maximum mean (max_mean) and median (max_med) response estimates (left and right, respectively). The minimum observed response (min_resp) is depicted with a blue circle and the maximum observed response value (max_resp) with a red circle (left lower). Finally, the minimum (min_logc) and maximum max_logc log10-scale concentrations are depicted, respectively, with blue and red vertical dashed lines (right lower).

- Concentration-Response Modeling Details

After summary values are obtained for each concentration-response series, all parametric models in tcplFit2 are fit to each series. Available model details are provided below:

Most models in tcplfit2 assume the background response is zero and the absolute response (or initial response) is increasing. In other words, these models fit a monotonic curve in either direction. The polynomial 2 (poly2) model is an exception with two parameterization options. The biphasic parameterization is what is used in tcpl . A biphasic poly2 model fits responses that are increasing first and then decreasing, and vice versa (assuming the background response is zero). If biphasic responses are not reasonable, data can be fit using the monotonic-only parameterization in a standalone application of tcplfit2_core with the parameter biphasic=FALSE assigned. This argument is not available in tcpl. All data is fit bidirectionally then responses in unintended direction may be indicated with negative hit calls if “overwrite” MC5 methods are applied.

Upon completion of model fitting, each model gets a success designation: 1 if the model optimization converges, 0 if the optimization fails, and NA if ‘nofit’ was set to TRUE within tcplFit2::tcplfit2_core function. Similarly, if the Hessian matrix was successfully inverted then 1 indicates a successful covariance calculation (cov); otherwise 0 is returned. Finally, in cases where ‘nofit’ was set to TRUE (within tcplFit2::tcplfit2_core ) or the model fit failed the Akaike information criterion (aic), root mean squared error (rme), model estimated responses (modl), model parameters (parameters), and the standard deviation of model parameters (parameter sds) are set to NA. A complete list of model output parameters is provided below:

Maximum likelihood estimation is utilized in the model fitting algorithm to estimate model parameters for all models. Even though tcplfit2 allows the maximum likelihood estimation to assume the error follows a normal or Student’s t-distribution, tcpl assumes the error always follows a t-distribution with four degrees of freedom. Heavier (i.e., wider) tails in the t-distribution diminish the influence of outlier values, and produce more robust estimates than the more commonly used normal distribution. Robust model fitting removes the need to eliminate potential outliers prior to fitting.

Let t(z, ν) be the Student’s t-distribution with ν degrees of freedom, y_i be the observed response at the i^th observation, and μ_i be the estimated response at the i^th observation. We calculate z_i as:

$$ z_{i} = \frac{y_{i} - \mu_{i}}{exp(\sigma)}, $$

where σ is the scale term. Then the log-likelihood is

$$ \sum_{i=1}^{n} [\ln\left(t(z_{i}, 4)\right) - \sigma]\mathrm{,} $$

where n is the number of observations.

The following plots provide simulated concentration-response curves to illustrate the general curve shapes captured by tcplFit2 models. When fitting ‘real-world’ experimental data, the resulting curve shapes will minimize the error between the observed data and the concentration-response curve. Thus, the shape for each model fit may or may not reflect what is illustrated below:

## Example Data ##
# example fit concentration series
ex_conc <- seq(0.03, 100, length.out = 100)

## Obtain the Continuous Fit of Level 4 Model Estimates ##
fits <- data.frame(
  # log-scale concentrations
  logc = log10(ex_conc),
  # parametric model fits from `tcplfit2`
  constant = tcplfit2::cnst(ps = c(er = 0.1), ex_conc),
  poly1 = tcplfit2::poly1(ps = c(a = 3.5, er = 0.1),x = ex_conc),
  poly2.mono.only = tcplfit2::poly2(ps = c(a = 0.13, b = 2, er = 0.1), x = ex_conc),
  poly2.biphasic = tcplfit2::poly2bmds(ps = c(b1 = 14, b2 = -0.1, er = 0.1), x = ex_conc),
  power = tcplfit2::pow(ps = c(a = 1.23, p = 1.45, er = 0.1), x = ex_conc),
  hill = tcplfit2::hillfn(ps = c(tp = 750, ga = 5, p = 1.76, er = 0.1), x = ex_conc),
  gnls = tcplfit2::gnls(ps = c(tp = 750, ga = 15, p = 1.45, la = 50, q = 1.34, er = 0.1),
                        x = ex_conc),
  exp2 = tcplfit2::exp2(ps = c(a = 0.45, b = 13.5, er = 0.1), x = ex_conc),
  exp3 = tcplfit2::exp3(ps = c(a = 1.67, b = 12.5, p = 0.87, er = 0.1), x = ex_conc),
  exp4 = tcplfit2::exp4(ps = c(tp = 895, ga = 15, er = 0.1), x = ex_conc),
  exp5 = tcplfit2::exp5(ps = c(tp = 793, ga = 6.25, p = 1.25, er = 0.1), x = ex_conc)
) %>% 
  reshape2::melt(data = .,measure.vars = c(
    "constant",
    "poly1","poly2.mono.only","poly2.biphasic","power",
    "hill","gnls","exp2","exp3","exp4","exp5"
  ))

## Updated Colors ##
fit_cols <-
  # choose 10 distinct colors
  viridis::magma(n = 11, direction = 1) %>% 
  # darken the original colors to make them more visible
  colorspace::darken(., amount = 0.2)

## Plot ##
fits %>%
  ggplot() +
  geom_line(aes(x = logc, y = value, lty = variable, colour = variable)) +
  facet_wrap(facets = "variable") +
  theme_bw() +
  labs(lty = "Models", colour = "Models") +
  scale_colour_manual(values = fit_cols) +
  ggtitle("General Shape of Models Included in `tcplfit2`") +
  xlab(expression(paste(log[10], "(Concentration) ", mu, "M"))) +
  ylab("Response")

This figure contains simulated concentration-response curves to illustrate the general underlying curve shape covered by each of the models included in the tcplfit2 package and used on the back-end of the level 4 data processing in tcpl. Each sub-plot in the figure corresponds to a single parametric model included in the model fitting process and has a corresponding color and line type to accompany it. All sub-plots are plotted such that the x-axis represents the log-transformed concentration (base = 10) and the y-axis represents the response values.

- Methods Assignment

For demonstrating the assignment of MC4 methods, method ‘1’ will be specified to calculate that the bmad and estimate one standard deviation of baseline using the two lowest concentration groups of treatment wells (well type, or wllt, equal to “t”).

Two examples of MC4 method assignments are provided below: (1) for a single assay endpoint and (2) all assay endpoints in containing string using tcplGetAeid . The assignment in the second approach can adapted for any subset of aeids.

## Methods Assignment #1 ##
# Assign the MC4 processing methods to aeid 80
tcplMthdAssign(
    lvl = 4, # processing level
    id = 80, # assay endpoint ID(s) to assign method(s)
    mthd_id = c(1), # method(s) to be assigned
    ordr = 1, # order the method(s) should be applied
    type = "mc") # the data/processing type

## Methods Assignment #2 ##
# Obtain the 'aeid' Values for all endpoints containing "ATG" string.
# "ATG" is the abbreviated assay source name of Attagene.
atg.aeid <- tcpl::tcplGetAeid(name = "ATG") 
## Assign the MC4 processing methods for subset of aeids 
tcpl::tcplMthdAssign( lvl = 4,
    id = atg.aeid[, aeid],
    mthd_id = c(1),  ordr = 1, type = "mc")

With the methods assigned, the MC4 processing can be completed for the desired set of aeids. After MC4 is processed, the user can load the model fit information from database.

A subset of MC4 data is available within the mc_vignette object.

The level 4 data includes fields for each of the ten model fits as well as the ID fields, as defined here. Model fit information are prefaced by the model abbreviations (e.g. cnst, hill, gnls, poly1, etc.). The fields ending in success indicate the convergence status of the model, where 1 means the model converged, 0 otherwise. NA values indicate the fitting algorithm did not attempt to fit the model. Smoothed model fits of the concentration-response data from the MC4 data object are displayed below:

# Load the example data from the `tcpl` package.
data(mc_vignette, package = 'tcpl')
# Allocate the level 3 example data to `mc3`.
mc3_example <- mc_vignette[['mc3']]
# level 3 does not store logc anymore, create it for plotting purposes
mc3_example[, logc := log10(conc)]
# Obtain the MC4 example data.
mc4_example <- mc_vignette[["mc4"]]

## Create a Sequence of Concentration Values within Observed Range ##
X <- seq(
 mc4_example[which(mc4_example[, spid] == "01504209"), conc_min],
 mc4_example[which(mc4_example[, spid] == "01504209"), conc_max],
 length.out = 100
)
## Obtain the Continuous Fit of Level 4 Model Estimates ##
# Apply each model fit to continous concentration values (X) and estimated
# parameters from 'tcplfit2'.
estDR <- mc4_example %>% 
  dplyr::filter(spid == "01504209") %>% 
  dplyr::reframe(
    cnst  = tcplfit2::cnst(.[, c(cnst_er)], x = X),
    poly1 = tcplfit2::poly1(.[, c(poly1_a, poly1_er)], x =  X),
    poly2 = tcplfit2::poly2(.[, c(poly2_a, poly2_b, poly2_er)], x = X),
    power = tcplfit2::pow(.[, c(pow_a, pow_p, pow_er)], x = X),
    hill  = tcplfit2::hillfn(.[, c(hill_tp, hill_ga, hill_p)], x = X),
    gnls  = tcplfit2::gnls(.[, c(gnls_tp, gnls_ga, gnls_p, gnls_la, gnls_q, gnls_er)], x = X),
    exp2  = tcplfit2::exp2(.[,c(exp2_a, exp2_b, exp2_er)], x = X),
    exp3  = tcplfit2::exp3(.[,c(exp3_a, exp3_b, exp3_p, exp3_er)], x = X),
    exp4  = tcplfit2::exp4(.[,c(exp4_tp, exp4_ga, exp4_er)], x = X),
    exp5  = tcplfit2::exp5(.[,c(exp5_tp, exp5_ga, exp5_p, exp5_er)], x = X) )
# Format data into a data.frame for ease of plotting.
estDR <- cbind.data.frame(X, estDR) %>%
  reshape2::melt(data = .,measure.vars = c(
    "cnst", "poly1", "poly2", "power", "hill", "gnls", "exp2", "exp3", "exp4", "exp5"))

## Updated Colors ##
fit_cols <-
  # choose 10 distinct colors
  viridis::magma(n = 10,direction = 1) %>% 
  # darken the original colors to make them more visible
  colorspace::darken(., amount = 0.2)

## Plot the Model Fits from Level 4 ##
mc3_example %>% 
  dplyr::filter(spid == "01504209") %>% 
  ggplot(.,aes(x = logc, y = resp))+
  geom_point(pch = 1, size = 2)+
  geom_line(data = estDR,
            aes(x = log10(X), y = value, colour = variable, lty = variable)) +
  labs(colour = "Models", lty = "Models") +
  scale_colour_manual(values = fit_cols) +
  xlab(expression(paste(log[10], "(Concentration) ", mu, "M"))) +
  ylab(expression(paste(log[2], "(Fold Induction)"))) +# )+
  ggtitle(
    label = paste("Level 4 Model Fits",
                  mc4_example[which(mc4_example[,spid] == "01504209"), dsstox_substance_id],
                  sep = "\n"),
    subtitle = paste("Assay Endpoint: ",
                  mc4_example[which(mc4_example[, spid] == "01504209"), aenm])) +
  theme_bw()

The plot depicts the observed concentration-response data with white circles, where the x-axis is base 10 log-transformed concentration values. All the ten model fits are displayed and distinguished by color and line-type.

How well the model is fitting the data (i.e. goodness of fit) can be approximated by the Akaike Information Criterion (AIC) and the root mean square error (RMSE or RME). For the AIC, let log(ℒ(θ̂, y)) be the log-likelihood of the model θ̂ given the observed values y, and K be the number of parameters in θ̂, then,

AIC = −2log (ℒ(θ̂, y)) + 2K.

The RMSE is given by

$$\mathrm{RMSE} = \sqrt{\frac{\sum_{i=1}^{N} (y_{i} - \mu_{i})^2}{N}}\mathrm{,}$$

where N is the number of observations, and μ_i and y_i are the estimated and observed values at the i^th observation, respectively.

> Level 5

Level 5 processing determines the winning model and activity for the concentration series, bins all of the concentration series into fitc categories, and calculates various potency estimates.

- Methods Assignment

The model with the lowest AIC value is selected as the winning model (modl ), and is used to determine the activity (or hit call) for the concentration series. If two models have equal AIC values, then the simpler model (i.e. model with fewer parameters) wins. Additionally, if the constant model describes the response best out of all the models fit then modl is reported as ‘none’. An example of side-by-side comparison of AIC values from each of the ten model fits for the example dataset mc_vignette is provided.

The summary values and estimated parameters from the winning model are stored in the respective mc5 and mc5_param tables. The activity of each concentration-response series is determined by calculating a continuous hit call that may be further binarized into active or inactive, depending on the level of stringency required by the user; herein, hitc < 0.9 are considered inactive. The efficacy cutoff value (coff) is defined as the maximum of all values given by the methods assigned at level 5. When two or more methods (i.e. cutoff values) are applied for processing, the largest cutoff value is always selected as the cutoff for the endpoint. In the event only one method is applied, then that will serve as the efficacy cutoff for the endpoint. Failing to assign a level 5 method will result in every concentration series being called active. For a complete list of level 5 methods, see tcplMthdList(lvl = 5) or ?MC5_Methods . See the Data Interpretation section for more details on hit calls and cutoff.

While the ToxCast pipeline supports bidirectional fitting, sometimes it is necessary to censor the hitc of curves fit in the biologically irrelevant direction. There are two methods for overwriting the hitc value, and if applied, these will overwrite the hitc value for any biologically irrelevant curve by flipping the hitc to a negative value.

The example we include in this vignette for demonstrating the assignment of level 5 methods specifies three different efficacy cutoff estimates for consideration. These efficacy cutoff estimates include 3 * bmad, log₂(1.2), and 5 * bmad, which correspond to mthd_id assignments 1, 3, and 5 respectively, and the largest of these three values will be selected as the cutoff for the endpoint. With the methods assigned, the data are ready for MC5 processing.

# Assign the MC5 Processing Methods to aeid 80
tcplMthdAssign(
  lvl = 5, # processing level
  id = 80, # assay endpoint ID(s) to assign method(s)
  mthd_id = c(1, 3, 5), # method(s) to be assigned
  ordr = 1:3, # order the method(s) should be assigned
  type = "mc") # the data/processing type

#MC5 processing for aeid 80
tcpl::tcplRun(
  id = atg.aeid[, aeid], # assay endpoint id to pipeline
  slvl = 5L, # level to start pipelining on
  elvl = 5L, # level to end pipelining on
  type = 'mc' # endpoint processing type - 'mc' = "multiple concentrations"
)

# Load MC 5 data for aeid 80
mc5 <- tcplLoadData(lvl = 5,
                    type = 'mc',
                    fld = 'aeid',
                    val = 80,
                    add.fld = TRUE)
# Prepare the data into a readable format.
mc5 <- tcplPrepOtpt(mc5)

A subset of MC5 data is available in the mc_vignette object, and includes fields for the best model fit, the potency estimates, other estimates from the best model fit, as well as the ID fields. The user can visualize the model fitting results using the tcplPlot functions. See the Data Retrieval with invitrodb and via API sections for more information.

# Allocate the level 5 data in `mc_vignette` to the `mc5` object.
mc5_example <- mc_vignette[["mc5"]]

For demonstrative purposes, an alternative visual of the model fits from MC4 and the best model as well as the potency estimates from MC5 data is produced below.{#mc5_plot}

## Obtain Data ##
# Load the example data from the `tcpl` package.
data(mc_vignette,package = 'tcpl')
# Allocate the level 3 example data to `mc3`.
mc3_example <- mc_vignette[['mc3']]
# level 3 does not store logc anymore, create it for plotting purposes
mc3_example[, logc := log10(conc)]
# Obtain the MC4 example data.
mc4_example <- mc_vignette[["mc4"]]
# Obtain the MC4 example data.
mc5_example <- mc_vignette[["mc5"]]
# First, we need to obtain the subset of data related to spid = "01504209",
# which is our example spid.
mc3_ss <- mc3_example %>% dplyr::filter(spid == "01504209") # Level 3 - conc-resp series
mc4_ss <- mc4_example %>% dplyr::filter(spid == "01504209") # Level 4 - model fits
mc5_ss <- mc5_example %>% dplyr::filter(spid == "01504209") # Level 5 - best fit & est.
# Next, we need to obtain the smooth curve estimate for the best model found
# in the Level 5 analyses of the `tcpl` pipeline.
# See Level 4 example above for how estDR is calculated.
estDR <- estDR %>%
  dplyr::mutate(., best_modl = ifelse(variable == mc5_ss[, modl],
                                     yes = "best model", no = NA))

## Generate a Base Concentration-Response Plot ##
basePlot <- mc3_ss %>% 
  # Observed Concentration-Response Data
  ggplot()+
  geom_point(aes(x = logc,y = resp),pch = 1,size = 2) +
  # Cutoff Band
  geom_rect(data = mc5_ss,
            aes(xmin = log10(conc_min), xmax = log10(conc_max), ymin = -coff, ymax = coff),
            alpha = 0.15, fill = "skyblue") +
  # Best Model Fit
  geom_line(data = dplyr::filter(estDR, variable == mc5_ss[,modl]),
            aes(x = log10(X), y = value,color = mc5_ss[,modl])) +
  scale_colour_manual(values = c("royalblue3"), aesthetics = "color") +
  # Other Model Fits
  geom_line(data = dplyr::filter(estDR,variable != mc5_ss[, modl]),
            aes(x = log10(X), y = value, lty = variable),
            alpha = 0.3, show.legend = TRUE) +
  # Legend Information
  labs(lty = "Other Models", color = "Best Fit") +
  # Titles and Labels
  xlab(expression(paste(log[10], "(Concentration) ", mu, "M"))) +
  ylab(expression(paste(log[2], "(Fold Induction)"))) +# )+
  ggtitle(
    label = paste("Level 5 Best Model Fit",
                  mc4_ss[which(mc4_ss[, spid] == "01504209"), dsstox_substance_id],
                  sep = "\n"),
    subtitle = paste("Assay Endpoint: ",
                     mc4_ss[which(mc4_ss[,spid] == "01504209"), aenm])) +
  # Background Plot Theme
  theme_bw()

## Potency Estimate Layers ##
# First, we need to obtain/assign colors for the potency estimates to be displayed.
potency_cols <-
  # choose 5 distinct colors
  viridis::plasma(n = 5, direction = -1) %>% 
  # darken the original colors to make them more visible
  colorspace::darken(., amount = 0.1)
  
## Compile the Full Level 5 Plot ##
linePlot <-
  # Start with the `basePlot` object.
  basePlot +
  # Next, add the various potency layers.
  # BMD
  geom_hline(
    data = mc5_ss,
    aes(yintercept = bmr),
    col = potency_cols[1]
  ) +
  geom_segment(
    data = mc5_ss,
    aes(x = log10(bmd), xend = log10(bmd), y = -0.5, yend = bmr),
    col = potency_cols[1]
  ) +
  geom_hline(
    data = mc5_ss,
    aes(yintercept = coff),
    col = potency_cols[2]
  ) +
  geom_segment(
    data = mc5_ss,
    aes(x = log10(acc), xend = log10(acc), y = -0.5, yend = coff),
    col = potency_cols[2]
  ) +
  geom_hline(
    data = mc5_ss,
    aes(yintercept = max_med * 0.5),
    col = potency_cols[3]
  ) +
  geom_segment(
    data = mc5_ss,
    aes(
      x = log10(ac50), xend = log10(ac50),
      y = -0.5, yend = max_med * 0.5
    ),
    col = potency_cols[3]
  ) +
  geom_hline(
    data = mc5_ss,
    aes(yintercept = max_med * 0.1),
    col = potency_cols[4]
  ) +
  geom_segment(
    data = mc5_ss,
    aes(
      x = log10(ac10), xend = log10(ac10),
      y = -0.5, yend = max_med * 0.1
    ),
    col = potency_cols[4]
  ) +
  geom_hline(
    data = mc5_ss,
    aes(yintercept = max_med * 0.05),
    col = potency_cols[5]
  ) +
  geom_segment(
    data = mc5_ss,
    aes(
      x = log10(ac5), xend = log10(ac5),
      y = -0.5, yend = max_med * 0.05
    ),
    col = potency_cols[5]
  )

# create data table for potency estimate points
mc5_points <- mc5_ss %>%
  select(bmd, acc, ac50, ac10, ac5) %>%
  tidyr::pivot_longer(everything(), names_to = "Potency Estimates") %>%
  mutate(x = log10(value)) %>%
  mutate(mc_color = potency_cols) %>%
  mutate(`Potency Estimates` = toupper(`Potency Estimates`))
yvals <- mc5_ss %>%
  select(bmr, coff, max_med) %>%
  tidyr::pivot_longer(everything()) %>%
  select(value) %>%
  mutate(reps = c(1, 1, 3)) %>%
  tidyr::uncount(reps) %>%
  mutate(y = value * c(1, 1, .5, .1, .05)) %>%
  select(y)
mc5_points <- mc5_points %>% cbind(yvals)

# add Potency Estimate Points and set colors
fullPlot <- linePlot + geom_point(
  data = mc5_points,
  aes(x = x, y = y, fill = `Potency Estimates`), shape = 21, cex = 2.5
) +
  scale_fill_manual(values = mc5_points %>% arrange(`Potency Estimates`) %>% pull(mc_color))

## Display the Compiled Plot ##
fullPlot

• Each of the concentration-response models fit in MC4 are included in the plot, where the blue curve indicates the best model fit for the observed data (white circles) and the rest are depicted by the gray curves. The light-blue shaded region represents the estimated efficacy cutoff (coff). The horizontal lines show the activity response levels from which potency estimates of interest are defined, and the vertical lines show the corresponding potency estimates. The black point shows the AC₅ (concentration producing 5% of the maximal response), the purple point shows the AC₁₀ (concentration producing 10% of the maximal response), the yellow point shows the BMD (benchmark dose), the orange point shows the ACC (concentration producing a response at the efficacy cutoff), and the pink point shows the AC₅₀ (concentration producing 50% of the maximal response).*

Additional information on derivations on potency estimates is found in Data Interpretation>Potency Estimates. After curve fitting, all concentration series are also assigned a fit category (fitc) based on similar characteristics and shape. See the Data Interpretation>Fit Category section for more details.

- Level 6

In addition to the continuous hitc and the fitc, cautionary flags on curve-fitting can provide context to interpret potential false positives (or negatives) in ToxCast data, enabling the user to decide the stringency with which to filter these targeted in vitro screening data. These flags are programmatically generated and indicate characteristics of a curve that need extra attention or potential anomalies in the curve or data. See the Data Interpretation>Flags section for more details.

- Level 7

For invitrodb v4.2 onward, a new mc7 table contains pre-generated AED values using several potency metrics from invitrodb and a subset of models from the High-throughput Toxicokinetics R package httk. AEDs are generated in a separate script using the httk R package. This is done separately due to the resource-intensive nature of running the Monte Carlo simulations to get estimates of plasma concentration for the median (50th %-ile) and most sensitive (95th %-ile) toxicokinetic individuals. Moreover, this is applied to both the 3-compartment steady state (3compartments) model and the physiologically-based toxicokinetic (pbtk) model for all chemicals included in invitrodb v4.2 (generation of the table as configured in the current code took 24h using 40 cores). See the Administered Equivalent Dose section.

Compiled Processing Examples

This section includes a practical applications for single- and multiple-concentration data, from methods assignment through data processing.

- SC Data

After SC0 is loaded into the database, SC1 and SC2 methods can be assigned. Once assigned, data can be processed.

## Methods Assignment
tcplMthdAssign(lvl = 1, id = 1:2, mthd_id = c(1, 11, 13), ordr = 1:3, type = "sc")
tcplMthdAssign(lvl = 2, id = 1, mthd_id = 3, type = "sc")
## SC0-2 Processing by acid
tcplRun(id = 1, type = "sc", slvl = 0, elvl = 2)

- MC Data

After MC0 is loaded into the database, methods can be assigned. Once assigned, data can be processed. This can be done in two ways: (A) from start to finish (i.e. Level 0 to 5) with the assay component ID (acid) or (B) Level 0 to 3 with the assay component ID (acid) and Level 4 to 5 with the assay endpoint ID (aeid). Option A may be helpful if data needs to be processed completely, whereas Option B or derivatives of Option B may be helpful when adjusting methods.

## Methods Assignment
# No MC1 methods needed
tcplMthdAssign(lvl = 2, id = 1, mthd_id = c(3,4,2), ordr = 1:3, type = "mc")
tcplMthdAssign(lvl = 3, id = 2, mthd_id = 1, ordr = 1, type = "mc")
tcplMthdAssign(lvl = 4, id = 2, mthd_id = 1, ordr = 1:2, type = "mc")
tcplMthdAssign(lvl = 5, id = 2, mthd_id = c(1,3,5), ordr = 1:3, type = "mc") 
tcplMthdAssign(lvl = 5, id = 2, mthd_id = c(1,3,5), ordr = 1:3, type = "mc") 

## Assign the Number of Processing Cores.
mycores <- 1 # If users do NOT want to leverage parallel computing. 
# Users that want to leverage parallel computing set to > 1, but less than the total number of cores
# (i.e. need at least 1 core open for overhead). If not provided, this will be assumed. 
# "parallel::detectCores()" can be run to understand the maximum allowed number of cores.

## Option A: MC0-5 Processing by acid
tcplRun(id = 80, type = "mc",  slvl = 0L,   elvl = 5L,  mc.cores = 20)

##Option B: MC0-3 Processing by acid, followed by MC4-6 by aeid
tcplRun(id = list$acid, type = "mc", slvl = 0L, elvl = 3L)
tcplRun(id = list$aeid, type = "mc", slvl = 4L, elvl = 5L)

Hit Calls

In tcpl v2, activity hit calls (hitc) were binary, where 0 was negative, 1 was positive, and −1 corresponded to concentration-response series that tcpl was “unable to fit” (e.g., <4 concentrations). In tcpl v3 onwards, the hitc is the product of three proportional weights, and the resulting continuous value is between 0 and 1, though the values are not normally distributed and tend to approach 0 or 1. Hitcalls that approach 1 indicate concentration-response series with biological activity in the measured response (i.e. ‘active’ hit). Hitc is typically binarized into active or inactive designations, depending on the level of stringency required by the user. For current ToxCast work, a hitc greater than or equal to 0.90 is labeled active, whereas anything less was considered inactive. This threshold of 0.90 was based on other tcplfit2 implementations with in vitro screening data (Nyffeler et al., 2023) and reflects the apparent bimodal nature of the hitc distribution, where a preponderance of the hitc fall between 0 and 0.1 and 0.9 and 1.0. Users may interpret the continuous hitc into active or inactive designations based on different thresholds. Further testing through implementation of this new functionality may reveal appropriate thresholds for different applications or assays. The “unable to fit” series now appear as model “none” with a hitc of 0 (inactive). Negative hitc in tcpl v3.2 correspond to curves that suggest biological activity in an unintended direction, as further described below.

Continuous hitc as defined in tcplfit2 R package is calculated as the product of three proportional weights representing the confidence that:

• p1: “the winning AIC value is less than that of the constant model.”
  • Determine whether the constant model – if allowed to win – is a better fit than the winning model – i.e., is the winning model essentially flat or not.The constant model may never be selected as the winning model, but if the constant model has the lowest AIC compared to other models, the calculated continuous hitc will be zero.
• p2: “at least one median response is greater than the cutoff.”
  • At least one dose group has a central tendency of the response values “outside” the cutoff band (consider bi-directional). Response is greater than cutoff in “+” direction and less than cutoff in “–” direction.
• p3: “the top of the fitted curve is above the cutoff”
  • Determine whether the predicted maximal response exceeds the cutoff, i.e. the response corresponding to the effect size of interest.

See Sheffield et al., 2021 for more information on tcplfit2.

Cutoff

The cutoff is a user-defined level of efficacy that corresponds to statistical and/or biological relevant change from baseline for each assay endpoint. All versions of tcpl provide methods for estimation of the baseline sampling variability, or noise around the assay controls, including calculation of the median absolute deviation over all response values given by wells that may represent baseline response (the BMAD), such as the neutral or vehicle control or the first two concentrations in the concentration series for all chemicals screened as defined by Level 4 methods. Users define mc5 methods depending on assay and data type, with some common cutoff thresholds used to establish a cutoff including 3 * BMAD, 20% percent change, or 1.2*log10 fold-change. Operationally in tcpl, the efficacy cutoff value (coff) is defined as the maximum of all values given by the methods assigned at level 5. When two or more methods (i.e. cutoff values) are applied for processing, the largest cutoff value is always selected as the cutoff for the endpoint. In the event only one method is applied, then that will serve as the efficacy cutoff for the endpoint. Failing to assign a level 5 method will result in every concentration series being called active. For a complete list of level 5 methods, see tcplMthdList(lvl = 5) or ?MC5_Methods.

Potency Estimates

Curve-fitting enables determination of various metrics of potency, i.e., concentrations at which some amount of in vitro bioactivity is expected to occur, as illustrated above. This includes Activity Concentrations at Specified Response and Benchmark Dose (BMD), which vary in the mathematical approach for computing these values, noting that logic for computation of the BMD is controlled in the R package tcplfit2.

- Activity Concentrations at Specified Response

An activity concentration is the estimated concentration inducing a specified level of response (activity). A common potency metric used from tcpl is the activity concentration at 50% of maximal activity, or AC50. The default baseline region is defined as ± 3 * BMAD⁠, and the ACB is the concentration at which the model first reaches a default of 3 * BMAD⁠, whereas ACC is defined as the concentration at which the model reaches the user-defined cutoff. Note that potency metrics such as AC50 are reported even if the hit call is inactive if a non-constant model could be fit to the data.

All versions of tcpl output the activity concentrations as described in Table 2.

- Benchmark Dose

A Benchmark Dose (BMD) is the activity concentration observed at the Benchmark Response (BMR) level. In the current implementation of tcpl and tcplfit2, BMR is only defined as 1.349 standard deviations of baseline response in the two lowest concentrations of treatment wells or neutral controls wells, as defined by Level 4 methods. tcpl uses the following definitions and assumptions for setting the BMR:

BMR is a change from the mean response at baseline (𝜇(𝑏)) by some multiple (𝑐) of the standard deviation of the baseline (𝑠𝑑(𝑏)).

Here, the baseline (𝑏) is defined as samples from the two lowest concentrations across chemicals within an assay endpoint and the 𝑐 = 1.349a (Yang et al., 2017).

A 90% confidence interval around the BMD, bounded by the benchmark dose lower bound (BMDL) and the benchmark dose upper bound (BMDU), is also computed and provided to reflect the uncertainty in the BMD estimate. The calculation of these confidence intervals will occasionally fail due to a singular matrix inverse, and in these cases, BMDU and BMDL will not be reported. This case occurs when the data are especially noisy and the confidence interval around the BMD approaches infinity. The winning model may return a bmd estimate that falls outside of the tested concentration range, so bounds are placed to censor the estimate values. The lower and upper bounds for bmd estimates are 0.1 * the lowest test concentration and 10 * the the highest test concentration, respectively. If the calculated bmd estimate is below or above the lower or the upper bounds, the value at the bound will be returned as the bounded bmd estimate instead.

Fit Categories

Fit Category Tree

A hierarchical fit category (fitc) decision tree is used to bin each fit as shown in Figure 2. Each fit falls into one leaf of the tree using the described logic with the final fitc indicated with gray boxes. Abbreviations are defined as: conc = concentration; hitc = hit call; |top| = absolute value of the modeled curve top; coff = cutoff; log_c(min) = minimum log₁₀ concentration tested; log_c(max) = maximum log₁₀ concentration tested; AC₅₀ = 50% activity concentration; AC₉₅ = 95% activity concentration.

After curve fitting, all concentration series are assigned a fit category (fitc) based on similar characteristics and shape. Logic is based on relative activity, efficacy, and potency comparisons as shown in Figure 5. For continuity purposes, fitc numbering has been conserved from past tcpl versions. Grouping all series into fitc enables quality control and can be useful in data cleaning applications, especially when considered with Level 6 flags. In invitrodb v3-3.5, a common filtering approach removed the least reproducible curve-fits, i.e. those with very low AC₅₀ (below the screened conc range) and low efficacy (within 1.2-fold of the cutoff) as well as 3+ flags. However, preliminary investigation into invitrodb v4.1-4.2 has suggested that removing curve fits with 4 or more flags, or possibly filtering based on specific flags in combination with fitc such as fitc 36, may be a more appropriate filtering approach due to changes in curve fitting and flags in invitrodb v4 and beyond. The stringency of filtering for flags should be explored in a fit-for-purpose way.

Fit category is largely based upon the relative efficacy and, in the case of actives, the location of the AC₅₀ and concentration at 95% activity (an estimate of maximum activity concentration, AC₉₅) compared to the tested concentration range. All concentration response curves are first split into active, inactive, or cannot determine. “Cannot determine” is indicative of exceptions that cannot be curve-fit, e.g. a concentration series with fewer than 4 concentrations. Active designations are determined for fitc based on whether the hitc surpasses the 0.90 threshold. For those series that are designated inactive with a hitc less than 0.90, fitc can be used to indicate to what extent the curve represents borderline inactivity via comparison of top modeled efficacy to the cutoff (i.e, the absolute value of the modeled top is less than 0.8 times the cutoff).

For active curves, efficacy, as represented by the modeled top, is compared to 1.2 times the cutoff (less than or equal to, or greater than), thereby differentiating curves that may represent borderline activity from moderate activity. Active curves also have potency metrics estimated, e.g., AC₅₀ and AC₉₅ values, that can be compared to the range of concentrations screened to indicate curves for which potency estimates are more quantitatively informative. Curves for which the AC₅₀ is less than or equal to the minimum concentration tested (fitc = 36, 40) may indicate AC₅₀ values that are less quantitatively informative than AC₅₀ values within the concentration range screened. When the AC₅₀ is greater than the minimum concentration tested but the AC₉₅ is greater than or equal to the maximum concentration tested (fitc = 38, 42), it is possible the maximum activity was not fully observed in the concentration range screened. Fitc for curves where the AC₅₀ and AC₉₅ are both within the concentration range screened (fitc = 37, 41) represent the most quantitatively informative AC₅₀ values.

Fitc 36 describes a curve that is of low efficacy and with a low AC₅₀, below the concentration range screened. These are more likely to be noise or less reproducible fits. Fitc 41 and 42 are the ideal fitc for reproducible curves, as demonstrated by these two fitc comprising the majority of positive (hitc > 0.9) curves in invitrodb v4.1. Fitc 40 indicates a curve with at least moderate efficacy, but an AC₅₀ below the concentration range screened. These chemicals may be positive or reference chemicals screened in the incorrect concentration window to observe their minimum activity. These curves may also represent high-confidence positives for which we have limited understanding of the slope of the concentration-response curve, and as such, the AC₅₀ may be associated with more uncertainty.

Flags (Level 6)

In addition to the continuous hitc and the fitc, cautionary flags on curve-fitting can provide context to interpret potential false positives (or negatives) in ToxCast data, enabling the user to decide the stringency with which to filter these targeted in vitro screening data. Cautionary flags on fitting were developed in previous versions of tcpl and have been stored at Level 6. These flags are programmatically generated and indicate characteristics of a curve that need extra attention or potential anomalies in the curve or data.

For example, a curve may be considered a single point hit with activity not at the highest concentration tested, but re-inspection of the flagged curve could indicate a potential false positive. Other flags may suggest borderline activity, overfitting, or cell viability assays that are fit with gain-loss as the winning model. It is important to note that flags have no effect on the hitc or potency estimates, but they may indicate that a curve requires further examination to aid in data interpretation. A full list of flags implemented and descriptions is presented below:

Representative Samples

Multiple samples of the same chemical (many spids per chid) may be tested for a given assay endpoint. The tcplSubsetChid function subsets multiple-concentration Level 5 (MC5) or single-concentration Level 2 (SC2) data to select the “best” single tested sample to represent that chemical for a given endpoint. The function uses a series of logic to select the “representative sample” in MC and SC. Representative sample designation for a tested chemical are stored in the MC5_chid and SC2_chid tables, where 1 indicates if ID (M5ID or S2ID, respectively) is a representative sample, else 0.

Variable Matrices

The tcplVarMat function creates chemical-by-endpoint matrices, combining summary information for single-concentration and/or multiple-concentration screening. Such a matrix can be useful for understanding which chemicals were tested, and what hit calls and potency values resulted from this screening. For example, one might ask, “Which chemicals were screened in either single-concentration or multi-concentration for some assay endpoint, and if MC screening was performed, what potency was estimated for any observed activity?”

For variable matrices from tcplVarMat , a representative sample is selected using the tcplSubsetChid function in cases where a chemical-by-endpoint was tested in multiple samples (performed separately for multiple-concentration or single-concentration screening). When multiple sample-assay series exist for one chemical, a single series is selected by the tcplSubsetChid function. See the Representative Samples section for more information.

1. “ac50” – The active concentration at 50% maximal response (AC₅₀) for the winning model.
   
2. “ac50_verbose” – The AC₅₀ for the winning model, with text inserted in place for some situations touched on below.
   
3. “acc” – The active concentration at user-defined cutoff (ACC ) for the winning model.
   
4. “acc_verbose” – The ACC for the winning model, with text inserted in place for some situations touched on below.
   
5. “mc_hitc” – The hit-call for the winning model in multiple-concentration (MC) screening.
   
6. “sc_hitc” – The hit-call in single concentration (SC) screening.
   

tcplVarMat produces matrices of combined sc-mc output. For the AC₅₀ and ACC matrices specifically, values are inserted in place to show complete views of what was tested and what the results. Further, “ac50_verbose” and “acc_verbose” replace such values with text. AC₅₀ and ACC are:

• Value as reported when the chemical is tested in MC and positive
• Set to 1e6 when the chemical is tested but negative in MC. In _verbose matrices, these are indicated as “MC neg”.
• Set to 1e7 when the chemical is not tested in MC but was screened in SC with a positive hitcall for the same aeid. In _verbose matrices, these are indicated as “SC pos, No MC”.
• Set to 1e8 when the chemical is not tested in MC but was screened in SC with a negative hitcall for the same aeid. In _verbose matrices, these are indicated as “SC neg, No MC”.
• Left as NULL if chemical is not tested in either MC or SC.

Both SC and MC data are required for tcplVarMat. As a result, the “API” driver is not currently supported for invitrodb v4.1 since it does not return SC data.

In addition, additional matrices can be defined by the ‘add.vars’ parameter. The ‘add.vars’ parameter will take any Level 4 or 5 field and create the respective matrix.

# create matrices with all chemicals and assays; the entire database
varmat <- tcplVarMat()
# create matrices using a subset of chemicals and/or aeids
aeids <- c(80)
dtxsid <- c("DTXSID80379721", "DTXSID10379991", "DTXSID7021106", "DTXSID1026081")
varmat <- tcplVarMat(aeid = aeids, dsstox_substance_id = dtxsid)
# create extra matrices by adding vars
varmat <- tcplVarMat(aeid = aeids, add.vars = c("m4id", "resp_max", "max_med"))

tcplVarMat returns a list of chemical by assay matrices (data.tables) where the rows are given by the dsstox_substance_id and corresponding chnm (chemical name) columns and the colnames are given by assay endpoint name (aenm). To export this list, save to a .xlsx. Each matrix will have its own sheet.

library(writexl)
write_xlsx(varmat, path = "varmat_output.xlsx")

Cytotoxicity Burst Distribution

Estimates of chemical concentrations that elicit cytotoxicity and/or cell stress have been informative for contextualizing bioactivity screening data in ToxCast by providing information on the likelihood that these data may be confounded by assay interference resulting from cytotoxicity and/or cell stress, particularly when a parallel or in-well estimate of cell viability is unavailable. As such, general estimates of the median and lower bound concentrations that might elicit cytotoxicity and/or cell stress in vitro have previously been calculated using the tcplCytoPt function, which considers activity across a suite of cell-based assays based on updates to previous work (Judson et al., 2016).

The burst threshold can be used to infer activity above or below an estimated threshold of generalized cell stress and/or cytotoxicity, where users can define the degree of difference between some observed bioactivity and estimates of cell stress and/or cytotoxicity needed in order to discern “selective” activity (i.e., bioactivity thought to occur independently of cell stress and/or cytotoxicity) vs. “non-selective” activity (i.e., bioactivity that appears to occur concomitantly with estimates of the concentration needed for cell stress and/or cytotoxicity).

These estimated concentration threshold values have been released in the “cytotox” table of invitrodb and are also provided on the CompTox Chemicals Dashboard (CCD) Bioactivity Summary Plot, as shown in the BPA example below.

tcplCytoPt function uses the assay endpoints listed in the burst_assay field of the “assay_component_endpoint” table, where 1 indicates including the assay endpoint in the calculation. The “burst” assay endpoints can be identified by running tcplLoadAeid(fld = “burst_assay”, val = 1) .

Assay Description Documents

Given ToxCast includes a heterogeneous set of assays across a diverse biological space, annotations in the database help flexibly aggregate and differentiate processed data whereas assay documentation aligned with international standardization efforts can make ToxCast data more useful and interpretable for use in decision-making. The OECD Guidance Document 211 (GD211) is a standard for comprehensive assay documentation describing non-guideline in vitro test methods and their interpretation. This template is intended to harmonize non-guideline, in vitro method descriptions to allow assessment of the relevance of the test method for biological responses of interest and the quality of the data produced. Unlike the assay element annotations which are often short in a standardized format or use a controlled term list, the assay_descriptions fields have no character limit for text. A compiled report of these assay description documents are available on the ToxCast Downloadable Data page.

Adminstered Equivalent Doses (Level 7)

The highest level assumption in the in vitro to in vivo extrapolation (IVIVE) approach employed here is that the in vitro bioactive concentration in a ToxCast assay endpoint is roughly equivalent to a human plasma concentration in vivo. For a review of IVIVE and httk models for it, please see: Breen et al, 2021

For invitrodb v4.2 onward, a new MC7 table contains pre-generated AED values using several potency metrics from invitrodb and a subset of models from the High-throughput Toxicokinetics R package httk . As implemented, this MC7 table provides users with pre-calculated estimates of the in vivo human administered dose (mg/kg/day) based on the in vitro bioactive concentrations as seen in ToxCast screening data.

Analytical QC and Applicability Domain

There is high value in understanding the outcomes of solubilization and chemical stability in the vehicle chosen to solubilize the chemical, i.e. a chemical’s applicability domain for in vitro screening. This informs what chemicals and samples should screened in future experiments. It also helps inform future structural models to understand which chemicals will be stable and detectable in solubilization, and further provide insight into possible degradation products that could be synthesized or purchased. Most critically, this information promotes understanding of uncertainty in estimates of initial experimental concentration of chemicals.

To establish a resource of applicability domain information at the substance and sample level, a retrospective analysis of the analytical QC data for the ToxCast/Tox21 chemical library was conducted. This involved reviewing legacy reports from gas chromatography-mass spectrometry (GCMS), liquid chromatography-mass spectrometry (LCMS), and Nuclear Magnetic Resonance (NMR) experiments. Additional Analytical QC, such as for the PFAS chemical library, and integration efforts are ongoing.

Retrieving Assay Element IDs

The tcplLoadAsid, tcplLoadAid, tcplLoadAcid, and tcplLoadAeid functions load relevant assay ids and names for the respective assay elements based on the user specified parameters.

# List all assay source IDs
tcplLoadAsid() 
# Create table of all assay endpoint ids (aeids) per assay source
aeids <- tcplLoadAeid(fld = "asid", # field to query on
                      val = 14, # value for each field
                      add.fld = c("aid", "anm", "acid", "acnm")) # additional fields to return

Retrieving Assay Annotations

Assay source, assay, assay component, and assay endpoint are registered via tcpl scripting into a collection of tables. The database structure takes the annotations and organizes them as attributes of the assay conductors, the assays (i.e., experiments), the assay components (i.e., raw readouts), or the assay endpoints (i.e., normalized component data) enabling aggregation and differentiation of the data generated through ToxCast and Tox21 (or other relevant partners). The annotations capture four types of information:

1. Identification information
2. Design information such as the technology, format, and objective aspects that decompress the assay’s innovations,
3. Target information, such as the target of technological measurement, biological intended target, and biological process, and
4. Analysis information about how the data were processed and analyzed.

# Select annotation and subset by ids or name, ex.
assay <- tcplQuery("SELECT * FROM invitrodb.assay where aid=1;")
component <- tcplQuery("SELECT * FROM invitrodb.assay_component;")
component <- subset(component, acid %in% source$acid)
endpoint <- tcplQuery("SELECT * FROM invitrodb.assay_component_endpoint;")
endpoint <- endpoint[grepl("ATG", endpoint$assay_component_endpoint_name),]

# Or select all annotations by joining multiple tables
annotations <- tcplQuery("SELECT * FROM invitrodb.assay
                  INNER JOIN invitrodb.assay_source on assay.asid=assay_source.asid
                  INNER JOIN invitrodb.assay_component on  assay_component.aid=assay.aid
                  INNER JOIN invitrodb.assay_component_endpoint on assay_component_endpoint.acid=assay_component.acid;")

Retrieving Chemical Information

The tcplLoadChem function returns all chemical information or can be filtered for user specified parameters, e.g. the chemical name (chnm) and chemical id (chid). The tcplLoadChemList function allows the user to subdivide the chemical IDs based on presence in different chemical lists. These chemical lists are curated by the US EPA in the Distributed Structure-Searchable Toxicity (DSSTox) database. Chemicals can belong to more than one chemical list, and will be listed as separate entries when loading chemical list information.

tcplLoadChem()
tcplLoadChemList(field = "chid", val = 1:2)

Retrieving Methods

The tcplMthdList function returns methods available for processing at a specified level (i.e. step in the tcpl pipeline). The user defined function in the following code chunk retrieves and outputs all available methods for both the SC and MC data levels.

# Create a function to list all available methods function (SC & MC).
method_list <- function() {
  # Single Concentration
  ## Level 1
  sc1 <- tcplMthdList(1, 'sc')
  sc1[, lvl := "sc1"]
  setnames(sc1, c("sc1_mthd", "sc1_mthd_id"), c("mthd", "mthd_id"))
  ## Level 2
  sc2 <- tcplMthdList(2, 'sc')
  sc2[, lvl := "sc2"]
  setnames(sc2, c("sc2_mthd", "sc2_mthd_id"), c("mthd", "mthd_id"))
  
  # Multiple Concentration
  ## Level 2
  mc2 <- tcplMthdList(2, 'mc')
  mc2[, lvl := "mc2"]
  setnames(mc2, c("mc2_mthd", "mc2_mthd_id"), c("mthd", "mthd_id"))
  ## Level 3
  mc3 <- tcplMthdList(3, 'mc')
  mc3[, lvl := "mc3"]
  setnames(mc3, c("mc3_mthd", "mc3_mthd_id"), c("mthd", "mthd_id"))
  ## Level 4
  mc4 <- tcplMthdList(4, 'mc')
  mc4[, lvl := "mc4"]
  setnames(mc4, c("mc4_mthd", "mc4_mthd_id"), c("mthd", "mthd_id"))
  ## Level 5
  mc5 <- tcplMthdList(5, 'mc')
  mc5[, lvl := "mc5"]
  setnames(mc5, c("mc5_mthd", "mc5_mthd_id"), c("mthd", "mthd_id"))
  # Compile the Output
  mthd.list <- rbind(sc1, sc2, mc2, mc3, mc4, mc5)
  mthd.list <- mthd.list[, c("lvl", "mthd_id", "mthd", "desc")]
  # Return the Results
  return(mthd.list)
}

# Run the 'method_list' functions and store output.
amthds <- method_list()
# Print the available methods list.
amthds

The tcplMthdLoad function returns the method assignments for specified id(s). Later sections provide more detailed examples for utilizing the tcplMthdLoad function for individuals ids.

Retrieving Data

The tcplQuery function allows a user to provide an SQL query to load data from the MySQL database into the R session. In the following chunk we provide an example, but any valid SQL query can replace the one provided. Please review Database Structure section to help construct these queries.

# Load sample table using a MySQL query.
samples <- tcplQuery("SELECT * FROM sample;")

The tcplLoadData function can be used to load the data from the MySQL database into the R session. Further, the tcplPrepOtpt function can be used in combination with tcplLoadData to add useful chemical and assay annotation information, mapped to the retrieved data.

When loading data, the user must indicate the applicable fields and ids for the corresponding data level of interest. Loading level 0 (SC0 and MC0), MC1, and MC2 data the assay component id (acid) will always be used. As described in Data Processing sections, SC1 and MC3 processing levels perform data normalization where assay component ids (acid) are converted to assay endpoint ids (aeid). Thus, the SC1 and MC3 data tables contain both acid and (aeid) ID’s. Data can be loaded using either id as long as it is properly specified. Loading SC2, MC4, 5, and 6 always use (aeid). Selected id(s) are based on the primary key within each table containing data.