MEDUSA for Identifying Death Regulatory Genes in Chemo-genetic Profiling Data

Summary

This protocol describes how to use the MEDUSA analytical method to quantify the death regulatory effect of each gene knockout. It includes instructions on determining experimental conditions that optimize sensitivity and a step-by-step tutorial on performing the analysis.

Abstract

Systematic screening of gain- or loss-of-function genetic perturbations can be used to characterize the genetic dependencies and mechanisms of regulation for essentially any cellular process of interest. These experiments typically involve profiling from a pool of single gene perturbations and how each genetic perturbation affects the relative cell fitness. When applied in the context of drug efficacy studies, often called chemo-genetic profiling, these methods should be effective at identifying drug mechanisms of action. Unfortunately, fitness-based chemo-genetic profiling studies are ineffective at identifying all components of a drug response. For instance, these studies generally fail to identify which genes regulate drug-induced cell death. Several issues contribute to obscuring death regulation in fitness-based screens, including the confounding effects of proliferation rate variation, variation in the drug-induced coordination between growth and death, and, in some cases, the inability to separate DNA from live and dead cells. MEDUSA is an analytical method for identifying death-regulatory genes in conventional chemo-genetic profiling data. It works by using computational simulations to estimate the growth and death rates that created an observed fitness profile rather than scoring fitness itself. Effective use of the method depends on optimal tittering of experimental conditions, including the drug dose, starting population size, and length of the assay. This manuscript will describe how to set up a chemo-genetic profiling study for MEDUSA-based analysis, and we will demonstrate how to use the method to quantify death rates in chemo-genetic profiling data.

Introduction

In chemo-genetic profiling, systematic genetic perturbations are used to understand the contribution of each gene to a given drug response^1,2,3. These experiments are valuable and can reveal important insights about drug responses, including the drug's binding target and mechanisms of drug influx/efflux. However, because these experiments are typically performed in a pooled manner that evaluates all genes simultaneously, generating mechanistic insights from chemo-genetic profiling data can be challenging.

The control mechanisms and genetic dependencies for drug-induced cell death tend to be challenging to resolve in chemo-genetic profiling data. There are several reasons for this problem, but many of these stem from the impacts of variation in cell proliferation⁴. For instance, because cells proliferate exponentially, a genetic perturbation's effect on the proliferation rate has a larger impact on the population size than changing cell death rates. Furthermore, because this biased sensitivity is exacerbated over time and because these experiments are typically performed over several weeks, most studies are optimized to be highly sensitive to proliferation defects and essentially insensitive to changes in drug-induced cell death. Other proliferation-related issues include varied coordination between proliferation and death (e.g., how fast is each clone growing while also dying, and does this vary between genetic perturbations) and the variations in proliferation rates for each clone in the absence of drug, which changes the expected number of cells that should/could have been recovered if the drug was not effective. The bottom line is that chemo-genetic profiling experiments generally score the impact of genetic perturbations using measurements that are proportional to the relative population sizes, comparing treated and untreated populations. Because population size is a product of both the cell growth and cell death rates, from the perspective of cell death, proliferation represents a confounding influence.

To remedy these issues, we created a Method for Evaluating Death Using a Simulation-assisted Approach (MEDUSA)⁵. MEDUSA works by interpreting the observed relative population size data through the lens of computational simulations to infer the combination of drug-induced growth and death rates that generated the observed drug response for each genetic clone. Prior data suggest that the method can accurately infer how genetic perturbations affect the drug-induced death rate, but the accuracy of this method depends on a detailed understanding of how cell proliferation and cell death are coordinated by a drug and how these rates vary over time^5,6. Additionally, MEDUSA-based inferences require a drug being tested at a dose that causes substantial cell death. Importantly, these drugging conditions create additional concerns about the starting population sizes and assay lengths, which should be carefully considered and optimized. In this protocol, we describe how to set up a chemo-genetic profiling study for MEDUSA-based analysis and provide a detailed use of this analytical method. The overall goal of MEDUSA is to determine how each gene deletion impacts drug-induced growth and death rates.

Protocol

1. Optimization of the drug dose to induce cell death

NOTE: The protocol below is for optimizing one combination of cell line, drug, and dose using the FLICK assay^7,8. The example data for this protocol describes the optimization and screening of 5 μM etoposide in U2OS cells, but the same optimization steps could be applied to any other desired combination of cell line and drug/dose. This protocol does not require the collection or analysis of dead cells. This protocol can be used for suspension cultures, provided that dead cells are removed. For suspension cultures, dead cells can be sorted/separated using a viability marker or a marker specific for cell death.

1. Before performing a FLICK assay, optimize cell number, Triton-X permeabilization, SYTOX concentration, and plate reader acquisition settings using the detailed protocol described in⁸.
2. Plate cells in black, clear-bottom 96-well plates. Typical seeding densities are 1500-5000 cells per well, and this number should be optimized per cell line. Incubate cells at 37 °C with 5% CO₂ and humidity for 16-24 h.
3. Prepare drug dilutions in media containing the concentration of SYTOX that was identified as optimal in step (1.1). Drugs are typically tested across a log dilution series containing 8-12 biologically or clinically relevant doses.
   NOTE: Certain drugs can emit fluorescence at the same wavelength as SYTOX, making quantification with this assay impossible. In such cases, utilizing SYTOX variants that emit fluorescence at different wavelengths might be beneficial. Additionally, this assay cannot assess the reaction to drugs that affect SYTOX fluorescence or prevent SYTOX from binding to DNA.
4. Lyse one plate (T0) using Triton-X at the concentration optimized in step (1.1). After cell lysis, use a plate reader to measure the SYTOX fluorescence and determine the total starting cell number based on the empirically established relationship between fluorescence and cell number (see FLICK protocol^7,8).
5. Measure SYTOX fluorescence in the experimental plates at T0 using a plate reader, with settings optimized in step (1.1).
6. Monitor the SYTOX signal over time for three days. To constrain the resulting kinetic data, use at least three time points per day, each 4 h apart.
   NOTE: This protocol evaluates cell death and growth over 3 days. This length of time is usually sufficient for death-inducing drugs; however, the assay can be shortened or lengthened to accommodate drugs with alternate death kinetics.
7. After the final timepoint, lyse the experimental plate(s) with Triton-X to determine the final population size of each condition.
8. Calculate Fractional Viability (FV) and fit endpoint data to dose curves, as described in the FLICK protocol^7,8.
9. Identify drug doses that induce roughly 50% cell death by evaluating the FV dose curves.
   NOTE: The identification of dead cells by SYTOX necessitates both plasma membrane rupture and nuclear envelope breakdown. These membranes need continuous upkeep, so SYTOX will identify dead cells regardless of the cell death mechanism induced by the drug. However, these processes do not happen with the same effectiveness for all forms of cell death, which could influence SYTOX's capacity to measure dead cells precisely. Additionally, different cell death mechanisms will produce cell corpses of varying stabilities. However, we have found that all forms of cell death tested so far result in SYTOX-bound DNA (either inside a permeabilized nucleus or in the cell culture media). This SYTOX signal is generally stable throughout a 72 h assay but can be monitored by evaluating the total cell number at the beginning and end of the assay.

2. Selection of assay length

1. Using the FLICK data collected in step 1, calculate each condition's Lethal Fraction (LF) over time as described in the FLICK protocol^7,8. Fit kinetic LF data to a lag-exponential death (LED) model, as described previously⁹.
2. Plot LF over time for each lethal dose identified in step 1. Select a dose and time point that produces ~50% LF.
   NOTE: A lethal fraction of ~50% enriches strong, death-specific signals in CRISPR screening data while leaving enough live cells behind to maintain library representation⁵. However, some drugs will be constrained by literature precedence or a narrow range of doses that induce a specific desired phenotype. If the dose is restricted, the assay endpoint should be shortened or extended to ensure the desired dose reaches ~50% LF. If 50% lethality is not possible, the size of the starting population may need to be increased.

3. Determining the starting population size

1. Select the desired drug dose to further optimize screening based on the above optimization.
2. Plate cells on 10 cm dishes in complete media. Each condition should be plated with technical triplicates and enough biological replicates so that cells can be counted every 12-24 h for the length of the screen, including the T = 0 time point.
   NOTE: Plated cell numbers may need to be optimized separately for untreated and treated cells, considering net population growth rate and cell size. Cells should generally be plated at densities that result in a confluency that does not compromise cell health at the assay endpoint (i.e., for many cell types < 80% confluent). Treated and untreated cells can also be replated during the screen, if necessary, as long as library representation is maintained.
3. Incubate cells overnight at 37 °C with 5% CO₂.
4. Add the drug at the desired concentration to the treatment plates. In parallel, harvest cells from 3 untreated plates and count the number of live cells using a hemocytometer. Stain the dead cells with trypan blue (or a comparable viability dye) to ensure that only live cells are counted. These data establish the number of cells at the T = 0 time point.
5. After 12-24 h, harvest cells from 3 untreated and 3 treated plates and count the number of live cells. Repeat the harvesting until the assay endpoint is reached. To fully understand the drug and cell death kinetics, it is recommended that live cells be monitored for 24-48 hours following the time point identified in step (2.2).
6. Visualize changes in viability by plotting population size (y-axis) over time (x-axis).
   NOTE: Based on the maximum observed cell number and the final population size, a rough estimate of the lethal fraction at the endpoint can be made. This should be compared to the fractional viability and lethal fraction data generated using the FLICK assay (step 1.9 and step 2.2) to ensure that the expected amount of cell death is achieved.
7. Identify the time point with the smallest population size. Depending on the amount of growth that occurs in the presence of the drug, this will be at either the start or end of the CRISPR screen.
8. Select a starting population size. Ensure the starting number of cells is optimized to maintain a minimum of 500x-1000x library coverage throughout the screen.
   NOTE: For example, the TKOv3 library has 70,948 sgRNAs. To represent this library at 500x coverage, no fewer than 35,474,000 cells per replicate per condition should be maintained throughout the assay. This bottleneck population size is often encountered at the beginning of the assay or following trypsinization and downsampling for untreated or proliferating conditions. In contrast, the bottleneck population size is often dictated by the population size at the assay endpoint for drug treatments that induce substantial lethality.
9. If desired, repeat the optimization above on the plates that will be used for the CRISPR screen (e.g., 15 cm plates) to confirm that drug efficacy scales correctly.

4. Performing CRISPR screen

NOTE: Cells can be screened using established CRISPR screening methods following optimization of drug dose and assay length. Established protocols, such as the TKO screening protocol from the Moffat lab^10,11, can be used to prepare the CRISPR screening library, produce viruses, perform the CRISPR screen, and prepare libraries for sequencing (Figure 1A-B). The protocol below describes steps specific to the MEDUSA analysis, which should be carried out on count-level data following sequencing.

1. Calculate experimentally observed fold changes for each sgRNA (Figure 1C) as described below.
   1. Generate a sgRNA-level counts table to prepare sequencing data for analysis. Trim and map sequencing libraries using standard pipelines, as described^5,12.
   2. Normalize sequencing library depth. Normalize libraries either manually or using built-in functions such as those in DESeq2, as in ^5,12. Perform normalization against all guides or by using the distribution of non-targeting guides.
   3. Randomly assign non-targeting sgRNAs to non-targeting genes. For example, a library with four sgRNAs per gene should have four non-targeting sgRNAs per non-targeting gene.
   4. For each comparison of interest (untreated/T0 and/or treated/untreated), calculate the log2(fold-change). It is recommended that this be calculated using a parametric fit in DESeq2, although a manual calculation of the fold change could also be performed. 
   5. Trim sgRNAs with a low number of counts. The correct trimming strategy will depend on the level of noise between replicates and how this varies over counts. These features will vary depending on the CRISPR library, and test multiple cutoffs in parallel for new data. Common cut-off strategies based on literature precedent include removing the lowest 5% of the guides or removing guides below some nominal threshold¹³.
2. Determine the impact of each single gene knockout on the growth rate as described below.
   1. Calculate the net population growth rate (NPG) of the untreated cells (i.e., the observed population doubling time, Figure 1D). This can be extracted from the FLICK data in step 1, live cell counts over time in step 3, or prior experimental measurements for the cell line of interest.
   2. MEDUSA requires three parameters to model the untreated condition: NPG: Net Population Growth rate, in doublings/h; T_start: Starting time point, in h. This is set to a default of 0; Tend_unt: Final time point for the untreated condition, in h.
   3. Simulate all possible perturbations to the NPG rate. Use a for loop to evaluate a range of possible NPG values. For each relative growth rate, calculate the log2(fold-change) that would be observed at the assay endpoint.
   4. For each experimental sgRNA in step 4.1.5, identify the simulated fold-change closest to the observed data. Assign the relative growth rate that resulted in the simulated fold-change to that sgRNA.
3. Characterize the coordination between growth and death as described below.
   NOTE: The parameters required for MEDUSA can be extracted from FLICK-style data when analyzed with the GRADE analysis method¹⁴. Data from (1) and (2) can be used for this, or an additional parameterization experiment can be performed.
   1. MEDUSA uses four parameters to characterize the coordination between growth and death in the presence of drug: GR_drug1: Drug-induced growth rate before death onset, in doublings/h; D_O: Death onset time, in h; GR_drug2: Drug-induced growth rate after death onset, in doublings/h; DR_drug: Drug-induced death rate, in units of Lethal Fraction/h. To parameterize D_O, use LF kinetics. Extract the death onset time from the LED fit.
4. Determine DR_drug using LF kinetics. Determine the average death rate by dividing the final LF by the time following death onset (in h).
5. Determine GR_drug1 from the live cell counting data in step 3. Before death onset, determine the drug-induced growth rate by fitting the data to a simple exponential equation.
   1. Determine GR_drug2 using a combination of data from step 3 and drug GRADE¹⁴. Calculate and plot the GRADE for the selected drug dose on a GRADE plot. If GRADE = 0, assume GR_drug2 to be 0; If GRADE > 0, use GRADE to determine the population's average growth rate. Based on the experimentally determined value for GR_drug1, determine the value for GR_drug2 that results in the observed average growth rate.
      NOTE: The average drug-induced growth rate can be inferred from a GRADE plot by computing the pairwise distance between the observed data point (i.e., pair of GR and FV data for a drug) and each point that makes up the GRADE space and identifying the closest point. For detailed instructions, see ref¹⁴. Drug-induced growth rates and death rates are not fixed values per drug and are specific to a drug at a given dose in a given cell context.
6. Generate simulated fold-changes for all possible combinations of growth rate and death rate as described below.
   1. Using the experimentally determined coordination of drug-induced growth rate and death rate, build a model that simulates the drug-induced population size over time (Figure 1D). This model can be built as a piece-wise function which combines two independent phases of the drug response (pre- and post-death onset time).
   2. For a full MEDUSA model, include these eight parameters: NPG: Net Population Growth rate, in doublings/h; GR_drug1: Drug-induced growth rate before death onset, in doublings/h; D_O: Death onset time, in h; GR_drug2: Drug-induced growth rate after death onset, in doublings/h; DR_drug: Drug-induced death rate, in units of Lethal Fraction/h; T_start: Starting time point, in h. This is set to a default of 0; Tend_unt: Final time point for the untreated condition, in h; Tend_tr: Final time point for the treated condition, in h.
   3. For the baseline model, simulate all possible growth and death rate perturbations. For each combination, calculate the log2(fold-change) that would be observed.
      NOTE: In MEDUSA, perturbations to the growth rate are assumed to be proportional across NPG, GR_drug1, and GR_drug2. For example, an 80% reduction in the NPG rate also results in an 80% reduction in both drug-induced growth rates.
7. Infer the death rate of each knockout (Figure 1D-E) as described below.
   1. For each sgRNA, triangulate the drug-induced death rate that would have produced the observed fold change. First, extract the relative growth rate determined in step (4.2.4).
   2. Determine rates for GR_drug1 and GR_drug2. Apply the relative growth rate from NPG to calculate a proportional growth rate for GR_drug1/GR_drug2.
   3. Using the constraints to NPG, GR_drug1, and GR_drug2, identify the simulated fold-change closest to the observed fold change for each sgRNA (treated/untreated). Use a combination of these observations to back-calculate the drug-induced death rate.
8. Calculate gene-level growth rates and death rates as described below.
   1. Determine gene-level rates for each gene in the library. Calculate the mean for each growth rate and death rate for all sgRNAs associated with a given gene (typically 4-10).
   2. To calculate empiric p-values, bootstrap the sgRNA-level data. Perform an FDR correction with the Benjamini-Hochberg procedure.
      NOTE: Complementary instructions, as well as sample data for analyzing CRISPR screening data using MEDUSA, are included on our GitHub repository (https://github.com/MJLee-Lab/MEDUSA). 

Results

Using this protocol, we explored the genetic dependencies of etoposide-induced death in U2OS cells. Etoposide is a commonly used DNA-damaging chemotherapy¹⁵. This chemo-genetic profiling experiment was performed using the GeCKOv2 library¹⁶ in Cas9 expressing U2OS cells⁵. In these types of experiments, gene function is typically inferred from the relative abundance of a given knockout clone when comparing drug-treated versus untreated populations (

Discussion

In this protocol, we describe the use of the MEDUSA analytical method for quantifying chemo-genetic profiling data to score how genetic perturbations alter drug-induced growth and death rates. Because the primary output of a chemo-genetic profiling experiment is related to population size, not rates, the underlying rates are inferred using simulations. Thus, the critical steps in the method are the accurate experimental parameterization of drug response kinetics, which includes the untreated cell proliferation rate, the ...

Materials

Name	Company	Catalog Number	Comments
100mm Fisherbrand TC Dishes	Fisher Scientific	FB012924	For growing cells and optimizing population size
150mm Fisherbrand TC Dishes	Fisher Scientific	FB012925	For growing cells and optimizing population size
DESeq2	R Package	v1.44.0	For calculating fold-change
DMEM	Corning	10017CV	For seeding and drugging cells
DMSO	Fisher Scientific	MT-25950CQC	For seeding and drugging cells
Fisherbrand 96-Well, Cell Culture-Treated, U-Shaped-Bottom Microplate	Fisher Scientific	FB012932	For seeding and drugging cells (pin plate)
Greiner Bio-One CELLSTAR μClear 96-well, Cell Culture-Treated, Flat-Bottom Microplate	Greiner	655090	For seeding and drugging cells
MATLAB	Mathworks	R2024a	For performing FLICK, GRADE, and MEDUSA
Microplate fluorescence reader	Tecan	Spark	For dead cell measurements
Sytox Green	Thermo Fisher Scientific	S7020	For dead cell measurements
TKOV3 CRISPR library	Addgene	125517	For performing CRISPR screen
Triton-X 100	Thermo Fisher Scientific	J66624-AP	For permeabilizing cells
Trypan blue	Corning	MT25900CI	For measure live/dead cells