Framework for in vivo T cell screens

CRISPR-based loss-of-function screens have enabled the simultaneous assessment of thousands of genes to discover new biological mechanisms and therapeutic targets in cancer cells (Shalem et al., 2014; Wang et al., 2014) and immune cells, including dendritic cells and T cells (Henriksson et al., 2019; Parnas et al., 2015; Shifrut et al., 2018). Thus far, most screens have been performed in vitro because, in that setting, it is more technically feasible to work with large numbers of cells and thereby perform large-scale screens (Doench, 2018). However, in vitro screens are limited in their ability to model multicellular interactions, complex differentiation states, and the metabolic milieu characteristic of in vivo immune responses (Chow and Chen, 2018; Miller et al., 2017). For example, T cell responses in vivo involve dendritic cell–mediated priming of naive T cells in the antigen-draining lymph node (LN), T cell differentiation into heterogeneous populations, trafficking to the target organ, and response to the inciting stimulus (e.g., pathogen or tumor) (Fowell and Kim, 2021). In vivo T cell screens in murine models have enabled interrogation of the regulation of complex phenotypes and identified genes important for T cell responses to tumors (Dong et al., 2019; Wei et al., 2019; Zhou et al., 2023) and chronic viral infection (Chen et al., 2021; LaFleur et al., 2019b), T cell differentiation into effector versus memory T cells (Huang et al., 2021), T follicular helper cell versus Type 1 helper cell differentiation (Fu et al., 2021), and inhibition of antigen-experienced CD4⁺ T cell proliferation (Sutra Del Galy et al., 2021).

Robust performance of in vivo screens requires adequate guide RNA (gRNA) representation—the presence of a sufficient proportion of the gRNA library and a sufficient number of cells per gRNA—to ensure a suitable sample size and replicability for downstream analyses (Dong et al., 2022). gRNA representation is, in turn, determined by the screen design parameters, including the size of the gRNA library and the number of replicate animals assessed and pooled (Doench, 2018). In addition, in vivo T cell screens typically involve adoptive transfer of ex vivo modified T cells into recipient animals followed by monitoring of T cell responses to an antigenic challenge. Accordingly, the number of T cells that can be adoptively transferred to new recipients before unduly perturbing the disease model as well as the engraftment rate of these cells into the recipients directly influence screen robustness. The majority of published in vivo T cell screens, including our own, did not describe prospective optimizations of screen design parameters (e.g., engraftment rate, mouse number, or library size) nor a retrospective analysis of the performance of the screen relative to the screen coverage (number of cells per gRNA), mouse number, or library size (Dong et al., 2019; Fu et al., 2021; Sutra Del Galy et al., 2021; Huang et al., 2021; LaFleur et al., 2019b; Wei et al., 2019; Zhou et al., 2023). Given the general lack of systematic experimental validation of in vivo screen design parameters, there is a need for guidelines to enable robust screens by providing methods for determining appropriate library size, cell transfer numbers, and numbers of mice that can be tailored to a broad spectrum of models.

One key readout of screens is the extent of enrichment or depletion resulting from a given gene knockout (KO); when applied to T cells, understanding the clonal nature of this expansion is often desired. Most in vivo T cell screens typically utilize TCR-transgenic T cells to control for antigen specificity, precluding use of the TCR to track clonal populations in vivo for understanding complex lineage relationships of T cells (Liu et al., 2022; Yost et al., 2019). This limitation has been circumvented in non-screen contexts through the use of sequenceable oligonucleotide barcodes (Gearty et al., 2022; Schepers et al., 2008), also known as unique molecular identifiers (UMIs) (Kivioja et al., 2011), that can be used to mark individual cells and their subsequent progeny. However, these powerful approaches have not yet been utilized in screening contexts in vivo to understand the impact of gene KOs on T cell clonal relationships. In addition, the phenotype of gRNA-containing cells within a population of cells is an aggregate of several disparate phenotypes: gene KO (intended effect), WT (unintended consequence of <100% editing efficiency), two gene KOs in one cell (unintended consequence of multiple integrants), and proliferative heterogeneity. UMIs can deconvolute the multiple possible phenotypes of gRNA-containing cells within a population of cells by discriminating between individual cell clones (Michlits et al., 2017; Schmierer et al., 2017; Zhu et al., 2019). While UMIs have been used to decrease the false-positive and false-negative rates of genome-wide in vitro CRISPR screens (Michlits et al., 2017; Schmierer et al., 2017; Xu et al., 2021; Zhu et al., 2019), UMIs have not been used to improve hit calling—the process of selecting genes for follow-up based on predefined statistical or effect size benchmarks—for in vivo CRISPR screens in T cells.

Here, we describe the Framework for In vivo T cell Screens (FITS), a framework for designing and analyzing in vivo T cell screens, that encompasses the following steps: (1) optimize cell transfer dose for a given cell type and model, (2) determine cell engraftment following cell transfer, (3) computationally model the optimal number of genes and mice needed for robust screening, (4) perform a targeted screen with UMIs, and (5) assess novel hits. As a proof-of-concept, we determined the optimal cell transfer dose and engraftment rate for CD8⁺ T cells in the tumor-draining LN (tdLN) and tumor using the B16-OVA (B16 melanoma tumor cells engineered to express the antigen ovalbumin) tumor model. Using this T cell dose and engraftment rate, we designed a 100-gRNA in vivo screen with 96 UMI-embedded tracrRNAs to assess genes that impact T cell proliferation in the tdLN and tumor as well as migration between the tdLN and tumor. This screen revealed site-specific regulators of T cell responses in tumors and tdLNs as well as putative regulators of T cell migration. Overall, FITS provides a template for designing and analyzing T cell CRISPR screens in vivo, which can be extended to other cell types and selection pressures.

We previously developed a cotransfer assay for assessing the cell-intrinsic effects of a given gene KO, compared with a control population, in naive CD8⁺ T cells responding to a tumor (LaFleur et al., 2019a). To assess CD8⁺ T cell–intrinsic regulators in tumors in a higher-throughput system, we developed an in vivo screening approach, optimizing three screen design parameters (Fig. 1 a): (1) the number of cells that can be adoptively transferred, (2) the number of genes that can be screened, and (3) the number of mice required to screen a given set of genes with a particular statistical power. To facilitate comparisons among distinct T cell clones, we restricted the T cell repertoire by using OT-1 T cell receptor transgenic mice (Hogquist et al., 1994) that recognize the peptide SIINFEKL from OVA protein expressed by B16-OVA tumor cells. We reasoned that adoptive transfer of supraphysiologic numbers of naive OT-1 CD8⁺ T cells would enable increased engraftment of T cells and therefore increase the number of gRNAs recovered in the tdLN and tumor. Transfer of too many naive OT-1 CD8⁺ T cells, however, could lead to tumor clearance, which would preclude analyses of tumor-infiltrating T cells. To determine the optimal number of naive OT-1 CD8⁺ T cells that could be adoptively transferred, we transferred a range of 20,000–300,000 cells and monitored tumor growth kinetics. Although adoptively transferred T cells slowed the growth of B16-OVA tumors, all transfer doses led to tumor outgrowth, suggesting the highest dose (300,000) would be optimal for screening (Fig. 1 b).

We next asked how many distinct T cell clones engrafted into recipient animals following an adoptive transfer of 300,000 naive OT-1 CD8⁺ T cells. We assessed this using a barcoding experiment on lentivirally transduced CD8⁺ T cells. Prior to performing the barcoding experiment, we confirmed that, following lentiviral transduction and culture, the CD8⁺ T cells were still naive (Fig. S1 a). We next performed a barcoding experiment by transducing naive OT-1 Cas9-expressing CD8⁺ T cells with a lentiviral library containing 6,000 distinct gRNAs (Table S1) and a Thy1.1 reporter to identify transduced T cells (Fig. 1 c). We performed these experiments with intergenic-targeting gRNAs to simulate the impact of double-strand breaks on cell fitness (Haapaniemi et al., 2018), and potentially subsequent engraftment, without directly targeting genes. To accurately measure clonal engraftment, it is important that each adoptively transferred T cell has a nearly unique genetic barcode. To accomplish this, we increased the complexity of this library to 576,000 by incorporating 96 UMIs (Michlits et al., 2017; Schmierer et al., 2017; Zhu et al., 2019) per gRNA into the tetraloop of the tracrRNA, which does not impact gRNA function (Zhu et al., 2019). Following transduction, we isolated live Thy1.1⁺ OT-1 Cas9⁺ CD8⁺ T cells (Fig. S1 b) to remove Thy1.1^– cells that may compete with the Thy1.1⁺ cells during engraftment into the tdLN and tumor. These Thy1.1⁺ cells were (1) transferred to Thy1.2⁺ recipient mice that express Cas9 to prevent immunologic rejection (Dubrot et al., 2021; LaFleur et al., 2019a) or (2) assessed for input representation of the gRNA-UMI combinations by next-generation sequencing. The next day, B16-OVA tumor cells were implanted into the recipient mice. After 8 days of tumor growth (Fig. S1 c), Thy1.1⁺ cells were isolated from the tumors and tdLNs (Fig. 1 d and Fig. S1 d) and analyzed by next-generation sequencing to determine the number of distinct gRNA-UMI combinations that were present. We estimated that ∼80% of the adoptively transferred cells would be clonal, with each gRNA-UMI represented by a single cell at the time of adoptive transfer (Fig. 1 e), suggesting that we can estimate T cell engraftment by counting the number of unique gRNA-UMI combinations.

The gRNA-UMI reads from the barcoding experiment span a range of values (Table S2), so the approach for counting a given gRNA-UMI as present will significantly affect apparent gRNA-UMI recovery. To determine the optimal read count cutoff, we examined how gRNA-UMI recovery changed with read count cutoffs from 0 to 100 and found that gRNA-UMI recovery plateaued for a read count cutoff of >1 read count (Fig. S1 e), suggesting this threshold may be sufficient for enumerating a stable number of recovered gRNA-UMI combinations. In addition, by randomly discarding reads, we found that the gRNA-UMIs recovered with a >1 read count cutoff were stable across downsampled sequencing read-depth (Fig. S1 f). Using this cutoff, on average 5,000 and 8,000 unique gRNA-UMI pairs were recovered in the tdLN and tumor, respectively (Fig. 1 f). Based on a 300,000-cell transfer, this represents an engraftment rate of ∼2% (Fig. 1 g). To validate that select gRNAs or UMIs were not impacting gRNA-UMI recovery, we assessed whether any of the 6,000 gRNAs or 96 UMIs had a consistent effect on T cell abundance in the tdLN or tumor across the replicate animals (Fig. S2, a–d). We found that the vast majority of the library is inert and can be used to assess the presence of distinct T cells from input through engraftment. To test if engraftment varies with the number of adoptively transferred T cells, we performed a second barcoding experiment examining a range of cell transfer doses (from 6,000 to 300,000 corresponding to 1 cell per gRNA to 50 cells per gRNA) and three timepoints after tumor implantation (Fig. 1 h and Fig. S2, e–g). Moreover, we included additional UMIs—384 instead of 96 (Table S3)—to increase the complexity of the library and thereby increase the fraction of cells that would be clonal at the time of transfer (Fig. 1 i). We found that the number of engrafted cells was relatively stable in the tdLN across cell dosages and timepoints (Fig. 1 j). In contrast, the number of engrafted cells in the tumor increased at later timepoints after tumor inoculation and with increasing cell doses, albeit with diminishing returns (Fig. 1 k). Thus, to achieve maximal engraftment, we chose 300,000 as our optimized cell transfer dose. Taken together, we optimized the cell transfer dose and determined, for this dose, the number of adoptively transferred CD8⁺ T cells that can engraft in the B16-OVA tumor model.

In vivo T cell screens enable evaluation of genes that control tissue-specific T cell responses, which cannot be studied in vitro. We hypothesized that different genes may uniquely impact T cell proliferation in the LN, LN-to-tumor migration, and intratumoral T cell proliferation (Fig. 2 a). To design our screen, we first evaluated (1) the library size that could be recovered, (2) the number of mice required to recover a given library size, and (3) whether it is preferable to perform several screens with smaller libraries in lower numbers of mice or a single larger library across more mice, given a fixed number of available mice. To answer these questions, we incorporated a requirement of three gRNAs per gene as most screen libraries typically contain multiple unique gRNAs per gene to account for variability in gRNA on-target efficiency and to control for off-target effects (Hanna and Doench, 2020). In addition, we mandated that each gRNA be recovered with a coverage (number of cells per gRNA) of 5 as it is important to have sufficient coverage of each gRNA to determine the biological effects of the gRNAs. Lastly, we simulated pooling mice together, as pooling leads to more efficient processing of tissues and is a common practice for in vivo screens (Fu et al., 2021; Huang et al., 2021; Wei et al., 2019). To calculate gRNA recovery in pooled samples, we assessed the recovery of a given number of gRNAs for each gene, requiring that the gRNAs would be present in an overall pool of 10 mice. We found that pooling 10 mice together and mandating recovery of 15 gRNAs per gene (to model three gRNAs each with a coverage of five cells) led to a recovery of ∼2,756 genes in the tdLN or tumor with a 90% probability of recovery (Fig. S3 a). In these calculations, the library size had a greater impact on library recovery than the number of mice (Fig. S3 a). To determine the library size we could screen for tdLN-to-tumor comparisons, we modeled the engraftment rate for tdLN and tumor each as 2%, corresponding to a lower-bound overall engraftment rate of both sites of 0.04%, based on the conservative assumption that the 2% rate applies both to initial tdLN engraftment and tdLN-to-tumor migration. Using this rate, we found that ∼29 genes could be recovered in 10 mice with a 90% probability of recovery (Fig. 2 b). Consistent with the tdLN and tumor analyses, for tdLN-to-tumor comparisons, library size impacts library recovery more than the number of mice used in the screen. Taken together, for a given number of available mice, it is preferable to perform several smaller screens in smaller sets of mice than a single large screen in all of the available mice. Furthermore, the constraints and goals of a given screen will affect the screen design. Thus, screeners can modify the variables of cell transfer dosage, engraftment rate, number of gRNAs per gene, coverage, and mouse pooling in these power calculations to design their screens using our publicly available scripts: https://github.com/scyrusm/fits.

Based on these power calculations, to enable evaluation of gene KO enrichment or depletion in the tdLN versus input, tumor versus input, and tumor versus tdLN, we designed a 100-gRNA library (Table S4) targeting 30 genes with three gRNAs per gene, within three gene classes: T cell coinhibitory/costimulatory receptors, T cell exhaustion, and inflammation sensing. In addition, we included 10 negative control gRNAs (five non-targeting and five intergenic-targeting) to benchmark the effects of the targeting gRNAs. To track cell replicates in vivo, we also included 96 UMIs per gRNA for an overall complexity of 9,600 gRNA-UMI combinations (Fig. 2 c). We adoptively transferred 300,000 transduced Thy1.1⁺ T cells to recipient mice that were subsequently injected with B16-OVA tumor cells. Following 11 days of tumor growth, we pooled the 48 mice with palpable tumors into six sets of eight mice each (pools A–F) and isolated transferred T cells from tumors and tdLNs (Fig. S3, b and c). To assess the representation of gRNAs, we performed next-generation sequencing of the input, tdLN, and tumor T cell samples, which revealed successful recovery of the majority of all possible gRNA-UMI combinations and all gRNAs in the tdLN and tumor (Fig. S3, d and e; and Table S5). Moreover, recovery of the gRNA library in the tdLN, tumor, or in the tdLN and the tumor (Fig. S3, f and g; and Fig. 2 d) fit well within the range predicted by our theoretical recovery calculations, despite using a smaller number of mice per pool (8 instead of 10). These data confirm our selection of (1) the library size and (2) the number of pooled mice required for successful library recovery.

Beyond simply recovering gRNAs, showing a consistent biological effect for a given gRNA is imperative for interpreting screen results. Therefore, we examined the consistency of the changes in gRNA abundance in the tdLN versus input, tumor versus input, and tdLN versus tumor. We found that our replicate groups of mice were consistent across all three comparisons (Fig. 2 e), suggesting we were recovering consistent biological effects. Thus, we next examined whether any gRNAs were enriched or depleted in the tdLN or the tumor. We first compared the fold changes of gRNAs for the tdLN versus input, tumor versus input, and tdLN versus tumor (Fig. S3 h; and Fig. 2, f and g). Examination of the tdLN to input revealed an enrichment of Zc3h12a- and Pdcd1-targeting gRNAs in the tdLN, as expected (Oh et al., 2020; Peng et al., 2020; Wei et al., 2019), whereas KO of Runx3 and Cd226 led to a depletion of T cells in tdLNs. In the tumor, Zc3h12a, Ptpn2, and Pdcd1 KO T cells were enriched, consistent with the known functions of these genes (Wei et al., 2019; LaFleur et al., 2019a). In addition, KO of Runx3, Il2ra, and Tbx21 led to a depletion of T cells in tumors. To identify site-specific regulators, we compared the tumor to the tdLN and found an enrichment of T cells with KO of Ptpn2 (LaFleur et al., 2019a), but not Pdcd1 and Zc3h12a (Fig. 2 g). Previous studies have shown an enrichment of Pdcd1 and Zc3h12a KO T cells in tumors, but a direct comparison of the magnitude of enrichment relative to the observed enrichment in the tdLN was not performed (Dong et al., 2019; Wei et al., 2019). Moreover, T cells with KO of Tbx21, Runx3, and Il2ra were depleted in the tumor relative to the tdLN, as previously shown for Runx3 KOs in the tumor versus the spleen (Milner et al., 2017). These data suggest Tbx21 may be a putative positive regulator of T cell trafficking to tumors, given the enrichment of Tbx21 KO T cells in the tdLN and depletion in the tumor as well as the known role of Tbx21 in inducing CXCR3 expression on CD8⁺ T cells (Taqueti et al., 2006). Collectively, these analyses confirm that the screen design parameters were appropriate for robust recovery of gRNAs, concordance among mice, and recovery of expected biology from our positive and neutral controls.

Given that the screen performed well, we retrospectively asked how the screen design parameters impacted the consistency of the screen results. To assess this, we downsampled the cellular coverage (cells per gRNA) by randomly excluding the read counts for a given number of UMIs because we previously showed that the UMIs were inert (Fig. S2, c and d). We considered decreasing numbers of UMIs (96 to 1) and assessed whether the gRNA fold changes in replicate sets of mice remained consistent (Fig. 3 a). This analysis revealed that the consistency of replicates increased with increasing number of UMIs analyzed and reached an intraclass correlation coefficient (ICC2) of 0.6 with ∼8 UMIs for the tdLN versus input and ∼20 UMIs for the tumor versus input. Given our theoretical calculation of ∼5× coverage per gRNA-UMI (see Materials and methods for details), this suggests that after engraftment ∼40 cells per gRNA for the tdLN versus input and ∼100 cells per gRNA for the tumor versus input is sufficient to achieve consistent gRNA fold changes across replicate animals (Fig. 3 b). In contrast, ∼75 UMIs (or ∼375 cells per gRNA) were required for the tumor versus tdLN to achieve an ICC2 of 0.6 (Fig. 3 c). Taken together, these data suggest that the coverage required for robust in vivo screens substantially changes for comparison between a given organ and an input sample or comparison between organs (Fig. 3 d).

UMIs can be used to track cells at monoclonal or oligoclonal resolution in screens (Michlits et al., 2017; Schmierer et al., 2017; Zhu et al., 2019). Each molecularly distinguishable clone can serve as an independent, biological replicate of a particular KO, increasing the effective sample size and markedly increasing statistical power to detect hits. Given that we included UMIs in our gRNA vector, we hypothesized that analyses of UMI replicates could be incorporated into assessments of effect size to further discriminate among putative hits and thus improve prioritization of genes for validation. Therefore, we calculated the fold changes of each gRNA-UMI combination in the tumor compared with the input (Fig. 4 a). We found that three of the four positive controls were still enriched (Zc3h12a and Pdcd1) and depleted (Runx3), as expected. However, we noted that (1) the effect sizes were highly variable even among the positive controls, for which we recovered the majority of the UMIs (Fig. 4 b) and (2) the rank order of putative hits changed relative to the gRNA level analysis. Thus, we asked whether a rank-based statistical test could be used to incorporate UMIs into hit calling. Previous studies developed the algorithms MAGeCKiBAR and ZFC (Z-score fold change) to analyze UMI level data in genome-wide in vitro CRISPR screens to decrease false-negative and false-positive rates (Xu et al., 2021; Zhu et al., 2019). We reasoned that these algorithms, which were benchmarked on genome-wide screens where the majority of genes will not have an effect, may not be appropriate for analyses of our targeted screen, where library composition may have an outsized impact on the significance of putative hits (Hanna and Doench, 2020). Indeed, analyses of all of the genes, just the enriched genes, or just the depleted genes with MAGeCKiBAR and ZFC revealed substantial differences in the statistical significance of putative hits (Table S6).

To address the dependence of hit significance on library composition, we developed a new algorithm, FAUST, which tests the null hypothesis H₀ that the mean proliferative effect of gRNAs targeting a particular gene is equal to that of a set of non-targeting and intergenic control gRNAs (Fig. 4 c). In brief, FAUST calculates the ratio R_i of counts for a gRNA-UMI i in an output sample versus an input sample, computes a Mann-Whitney U statistic (by default, though FAUST allows custom parametric tests to be used) for R_i associated with experimental versus control targets, and derives a P value corresponding to the above null hypothesis via the Mann–Whitney U test. The application of the Mann–Whitney U test provides two advantages to the parametric tests used previously. First, gRNA-UMI count tables are characterized by strongly enriched outliers (Fig. S4 a and Table S5), leading to distributions that are poorly modeled by a negative binomial distribution (Fig. S4 b), which is an assumption of the MAGeCKiBAR (Zhu et al., 2019) and ZFC (Xu et al., 2021) algorithms. In contrast, the Mann–Whitney U test requires no parametric assumptions, is robust to these outliers, and is computationally efficient. Secondly, the application of the Mann–Whitney U test produces values that are agnostic to sample-level normalization. Therefore, sample-level differences in read-depth or engraftment (which affect experimental and control gRNA-UMIs proportionally) will not affect whether genes are detected as hits. Having computed P values in this manner for each replicate, we can aggregate replicate pools of mice in one of two ways. We can either test the union null hypothesis (that the null hypothesis H₀ is true for at least one of the replicates) or the intersection null hypothesis (that the null hypothesis H₀ is true for all of the replicates). We address the former—hereafter called the conservative approach—by taking the maximum P value across all replicates. We address the latter—hereafter called the liberal approach—by aggregating P values for all replicates by Fisher’s combined test. These aggregate p-values are then adjusted to account for multiple hypothesis testing by applying the Benjamini-Hochberg procedure (q-values indicating an overall false discovery rate).

We applied FAUST to the 9,600-gRNA-UMI targeted screen. We first assessed the common language effect size (CLES), the probability of recovering a given target gRNA or a control gRNA following a random sampling. Importantly, the CLES of R_i was consistent across biological replicates despite differences in read depth and cell recovery (Fig. 4, d–f). Comparison of gRNA representation in the tdLN to input (Fig. 4 d), using the conservative mode of FAUST, revealed 11 putative hits (q value <0.05): Zc3h12a, Pdcd1, Tbx21, Eomes, Ptpn2, Cblb, Ptpn1, Fli1, Tcf7, Cd226, and Runx3. Thus, we were able to recover our positive controls (Zc3h12a, Pdcd1, and Ptpn2) (Oh et al., 2020; Peng et al., 2020; Wei et al., 2019) as well as previously unknown regulators of tdLN accumulation (Tbx21 and Ptpn1). Moreover, FAUST identified 15 genetic targets that significantly regulated T cell proliferation in the tumor, including our positive controls (Zc3h12a, Ptpn2, Pdcd1, and Runx3) (Fig. 4 e). To assess putative regulators of migration or site-specific responses, we compared the tumor with the tdLN and identified 8 of 30 genetic targets that significantly regulated T cell proliferation in tumors compared with tdLN, with one positive and seven negative regulators (Fig. 4 f). We confirmed the prior fold change analysis implicating Tbx21 as a putative positive regulator of T cell trafficking to tumors (Fig. 2, f and g; and Fig. S3 h), given the enrichment of Tbx21 KO T cells in the tdLN and depletion in the tumor. The stringency of calling a putative hit depends on the number of genes that can feasibly be subsequently validated. Thus, we also employed the less-stringent liberal mode of FAUST (Fig. S4, c–e), which yielded many additional putative hits, particularly for the tumor versus tdLN comparison. Tcf7 KO T cells were enriched in the tumor compared with the tdLN (Fig. S4 e), concordant with a requirement for TCF1 in stem-like CD8⁺ T cells (Kurtulus et al., 2019), which are predominantly found in tdLNs (Connolly et al., 2021). To determine the applicability of FAUST for screens that do not incorporate UMIs, we analyzed the 100-gRNA screen at the gRNA level using FAUST or MAGeCK (Li et al., 2014) and found that FAUST identified more enriched and depleted hits (Fig. S4, f–i). Collectively, through benchmarking to internal controls, FAUST improves hit calling for targeted screens.

Beyond gRNA enrichment or depletion in a given organ, we hypothesized that, for individual gRNAs, UMIs could be used for tracking T cell responses primed from multiple LN sites that migrate to tumors. For tumors injected on the ventral flank, the ipsilateral axillary and inguinal LNs will drain tumor antigens (Fransen et al., 2018; Marzo et al., 1999), enabling evaluation of antigen-specific T cell responses. In addition, we observed expression of CD44 and PD-1 on OT-1 CD8⁺ T cells in the contralateral LNs following B16-OVA challenge (Fig. S5 a), suggesting (1) antigen drainage and T cell priming in these LNs or (2) recirculation of cells originating from the ipsilateral LNs. Thus, we repeated the 100-gRNA screen with 96 UMIs per gRNA examining T cell responses in the tumor and these four LNs (Fig. 5 a). Following 11 days of tumor growth (Fig. S5 b), we sorted transferred cells from the tumor and four LNs (Fig. S5 c) and examined gRNA-UMI recovery. We recovered the majority of the 9,600 possible gRNA-UMI combinations in the tumor and ipsilateral LNs (Fig. S5 d). However, in the contralateral inguinal and axillary LNs, we recovered 2,648 and 1,967 gRNA-UMIs, respectively, consistent with the lower number of transferred cells isolated (Fig. S5 c). We next utilized FAUST to examine the CLES in the tumor and LNs and found that overall the CLES for the 30 gene KOs in the library were most similar among LNs (Fig. 5 b and Table S7). Some gene KOs, such as Pdcd1 and Zc3h12a, were enriched in all of the LNs and the tumor, whereas other gene KOs, such as Ptpn1, Cblb, and Ptpn2, were more enriched in the tumor compared with the LNs. Moreover, Tbx21 and Eomes KOs were most enriched in the ipsilateral LNs compared with the contralateral LNs, which may be due to (1) the varied tumor antigen load and (2) a concentration gradient of soluble factors, such as cytokines, at proximal compared with distal sites. Collectively, this screen identified gene KOs that have effects in LNs broadly, tumor-proximal LNs, or tumors, suggesting site-specific regulation of T cell responses to tumors.

To understand the site-specific phenotypes of gene KO T cell clones further, we assessed the degree of sharing of T cell clones by tracking UMIs between multiple LN sites and the tumors using Morisita’s overlap index (MOI), a statistical measure used to compare the overlap among two samples (Fig. 5 c). Based on our screen parameters (300,000 cell transfer, pooling of nine recipient animals, a 2% engraftment rate, and 9,600 gRNA-UMI combinations), each gRNA-UMI should be represented by approximately five cells, suggesting we can examine responses at oligoclonal but not monoclonal resolution in this screen. Assessing the MOI of the LNs and tumors, we found that the ipsilateral inguinal and axillary LNs shared a large number of clones for the majority of gene KOs assessed. In contrast, the tumor shared a large number of clones with the ipsilateral inguinal and axillary LNs for a subset of gene KOs: Eomes, Ptpn2, Ptpn1, and Cblb. Different gene KOs showed considerable heterogeneity in MOI across different tissues (Fig. 5 c): Zc3h12a and Pdcd1 KOs showed high LN MOI across all four LNs, whereas Tbx21 KOs had high ipsilateral LN MOI. We next examined whether the magnitude of enrichment or depletion of a gene KO correlated with shared clonality across the periphery. We reasoned that the average MOI for all pairwise LN comparisons would serve as a surrogate for similarity of T cell clonal responses for a given gene KO in the periphery (Fig. S5 e). Comparison of tumor versus input CLES and LN MOI revealed a dichotomy in T cell responses where some gene KOs, such as Zc3h12a and Pdcd1, had a high LN MOI and a high tumor CLES, whereas other gene KOs, such as Cblb, Ptpn2, and Ptpn1, had a low LN MOI and a high tumor CLES (Fig. 5 d). Thus, the tumor versus input CLES and LN MOI were not correlated. In contrast, the CLES for all four LN versus input comparisons correlated with LN MOI (Fig. 5 e and Fig. S5, f–h). These data suggest that there may be several paths for a CD8⁺ T cell to become enriched in the tumor: (1) high global expansion of T cells in LNs and tumors, such as Zc3h12a KO T cells, or (2) high intratumoral expansion, such as Ptpn2 KO T cells. To assess this hypothesis, we analyzed the read counts of individual UMIs for gRNAs targeting Zc3h12a and Ptpn2 in the input, LN, and tumor samples. We found that the majority of Zc3h12a KO T cell clones were expanded in the LN and tumor, whereas Ptpn2 KO T cell clones were primarily expanded in the tumor (Fig. 5 f and Fig. S5 i). Thus, analyses of individual UMIs revealed expansion kinetics of clonal populations and gene KO-specific effects on expansion dynamics. Collectively, analyses of T cell clonality, using UMIs, suggest that there are at least two paths—global expansion and intratumoral expansion—that result in CD8⁺ T cell enrichment in tumors (Fig. 5 g).

CRISPR screens have been widely performed in vitro and have well-established screen design rules. However, optimal translation of these parameters to in vivo screens, many of which depend on adoptive transfer of a pool of modified cells, is complicated by limited cell engraftment and mouse-to-mouse heterogeneity, necessitating new tools to aid in screen design. Here, we developed FITS, a framework for experimentally determining in vivo screen parameters and analyzing screens, using CD8⁺ T cells as a model cell type (Fig. 6). We demonstrate the utility of FITS by designing and executing a screen of 9,600 gRNA-UMIs in CD8⁺ T cells to define regulators of T cell expansion in response to a tumor challenge through paired analysis of the tumor and multiple tdLNs. Using gRNA-incorporated UMIs as molecular barcodes, we identified site-specific regulators of antitumor T cell responses with LN-dominant effects of Pdcd1 and Zc3h12a gene KOs and tumor-dominant effects of Ptpn1 and Ptpn2 gene KOs. Thus, FITS provides a template for designing and analyzing robust in vivo T cell screens in additional immunologic contexts.

Given the reliance of most in vivo screens on adoptive transfer of ex vivo manipulated cells to a new host, the engraftment rate of these cells greatly impacts the number of modified cells that can be assessed. Here, we performed an in vivo barcoding experiment in the tdLN and tumor and recovered ∼2% of our transferred cells. This number is likely model and cell-context specific, highlighting the importance of initial barcoding experiments in our framework (Fig. 6). For example, the strength of selection pressure, such as tumor model immunogenicity and pMHC-TCR affinity, influences the signal-to-noise of a given screen (Pan et al., 2018; Zhou et al., 2014), supporting the need to characterize cell engraftment and gRNA recovery more broadly for each cell type and disease model. Moreover, following adoptive transfer of CD8⁺ T cells, it is unknown whether an absolute number or a fixed percentage of cells engraft. Here, we addressed this question by performing a barcoding experiment with a fixed library size and a titration of adoptive cell transfer dosages. We found that transfer of increasing numbers of cells only marginally increased engraftment numbers, which suggests that an absolute number of cells engraft. Limited by overall engraftment rates, in vivo screens often employ pooling of mice together to facilitate tissue processing and to improve recovery of gRNAs. However, whether it is optimal to perform multiple smaller screens in smaller sets of mice or a single larger screen in a larger set of mice is an open question. We modeled recovery across a range of library sizes and numbers of mice to determine the optimal number of gRNA-containing cells we could reproducibly assess. Interestingly, changing the number of mice did not impact the likelihood of recovering gRNAs as much as changing the number of gRNAs assessed. Thus, in our model, it is more favorable to perform multiple smaller screens as opposed to a single large screen. This has additional advantages as it better controls for mouse-to-mouse heterogeneity, which we found to significantly impact engraftment rates. Following engraftment, the coverage needed for robust in vivo screens is unknown. To assess this, we downsampled our screen data to determine the optimal coverage for recovering concordant gRNA fold changes and found that ∼40 cells per gRNA were required for tdLN versus input analyses, whereas ∼100 cells per gRNA were required for tumor versus input analyses. As expected, interorgan (e.g., tumor versus tdLN) analyses required greater coverage (∼375 cells per gRNA). We expect the optimal coverage to also be model specific; the analyses presented here provide a useful framework for empirical determination of this value in other contexts.

UMIs are a powerful tool for decreasing false-positive and false-negative rates in screens (Michlits et al., 2017; Schmierer et al., 2017; Xu et al., 2021; Zhu et al., 2019), but they have been primarily used in vitro where it is straightforward to control experimental conditions and screen coverage. In contrast, in vivo experiments are subject to bottlenecking, which may impact false-positive and false-negative rates, and ultimately restricts screen size. We reasoned that incorporation of UMIs, which can serve as distinct cellular replicates for each given gRNA, could provide the necessary statistical power to overcome the challenges of hit calling unique to targeted in vivo screens. Previous approaches to identify gene-level drivers or inhibitors of cell proliferation compared the effect of a gRNA targeting a particular gene to the effect of all other genes in the screen (Li et al., 2014; Xu et al., 2021; Zhu et al., 2019) and may normalize for sample-level differences in read depth via the median-ratio method. Both of these strategies, while appropriate for identifying dominant drivers of proliferation amidst a large panel of targets with no effect, cause a target’s apparent effect size and significance to be skewed by other targets in the screen, which would particularly bias analyses of curated gRNA libraries. We implemented a method for targeted screen analysis (FAUST), comparing target gRNAs to predefined negative control gRNAs, circumventing the effect of library composition and size on statistical assessment of putative hits. This approach incorporates UMIs for a given gRNA as explicit cellular replicates to increase statistical power, in contrast to other algorithms that implement UMIs to provide a variance adjustment for each gRNA (Michlits et al., 2017; Xu et al., 2021; Zhu et al., 2019). Using our approach, we were able to recover known and novel regulators of T cell abundance in tumors and tdLNs with remarkable consistency between mouse biological replicates. To examine the broad utility of FAUST beyond UMIs, we compared hit calling at the gRNA level with FAUST and MAGeCK. We found that FAUST outperforms MAGeCK for both enriched and depleted hits in our 100-gRNA screen. Collectively, FAUST fills a gap in the suite of screen analysis algorithms by enabling more powered analyses of targeted screens and improving incorporation of UMI-level data as cellular replicates.

In addition to improving statistical power and enabling comparisons across screens, incorporation of UMIs increases the depth of information that can be extracted from a single screen. Hypotheses can be generated about how a given gene KO impacts peripheral versus local responses. As a proof-of-concept, we utilized UMIs to track T cell responses in tumors and multiple LN sites to identify gene KOs that lead to global cell expansion (e.g., Zc3h12a) and those that lead to local cell expansion (e.g., Ptpn2). One limitation of our analyses using UMIs to track clonal responses is the number of UMIs utilized (96) did not enable strictly monoclonal analyses. We envision incorporating a more diverse set of UMIs, such as the 384 UMI library used in Fig. 1, into future screens will enable monoclonal analyses of T cells with distinct gene KOs. Furthermore, to increase the power of these approaches, future studies combining screens, UMIs, and single-cell RNA sequencing will enable screeners to track the effects of gene KOs on clonal relationships and cell states (Zhou et al., 2023). Taken together, UMI-embedded tracrRNAs are a powerful addition to in vivo screens by enabling robust hit calling and assessment of clonal responses.

In summary, FITS serves as a novel experimental and analytical framework for designing and analyzing in vivo T cell screens. Here, we provide a proof-of-principle of the steps and analyses a screener could perform to ensure robust recovery and hit calls from a targeted in vivo screen. Moreover, we introduce the use of UMIs for cellular replicate analyses as well as for generating hypotheses about how gene KOs impact cell clonality and site-specific responses. Lastly, although we focused on in vivo T cell screens here, we envision that this framework of (1) optimizing cell transfer dose, (2) determining cell engraftment following cell transfer, (3) computationally modeling the optimal number of genes and mice needed, (4) performing a targeted screen with UMIs, and (5) assessing novel hits could be used for in vivo screens in other disease contexts and immune cell types.

7–16-wk-old male and female Cas9-expressing mice were used as recipients of T cells for all experiments, and 5–29-wk-old male and female CD45.1, CD45.1/CD45.2, or CD45.2 OT-1 Cas9-expressing mice were used as spleen donors for barcoding and screen experiments. LoxP-STOP-loxP-Cas9 mice (B6J.129(B6N)-Gt(ROSA)26Sor^tm1(CAG-cas9*^,-EGFP)Fezh/J Jax# 026175) (Platt et al., 2014) were backcrossed 10+ generations to C57BL/6J mice (# 000664; Jax) and then bred to Zp3-Cre mice (C57BL/6-Tg[Zp3-cre]1Gwh/J Jax# 003651) to delete the loxP-STOP-LoxP in the female germline. We subsequently bred out the Zp3-Cre allele. This Rosa26-Cas9–expressing strain (LaFleur et al., 2019b) was then bred to CD45.1 mice (B6.SJL-Ptprc^aPepc^b/BoyJ Jax# 002014) for marking transferred cells. The CD45.1 Cas9-expressing strain was then bred to OT-1 mice (C57BL/6-Tg[TcraTcrb]1100Mjb/J Jax# 003831) (Hogquist et al., 1994) to restrict TCR responses to the OVA antigen SIINFEKL. Mice were genotyped using PCR or Transnetyx probes.

Age- and sex-matched animals were used for each experiment. Within a given experiment, Rosa26-Cas9–expressing mice were sourced from in-house breeding colonies and Jackson Laboratories (# 026179; Jax). Mice were cohoused to control for potential differences in their microbiomes. All experimental mice were housed in specific pathogen–free conditions with a 12-h light cycle and a 12-h dark cycle, 30–70% relative humidity, and 21.1°C ambient temperature. Mice were used in accordance with animal care guidelines from the Harvard Medical School Standing Committee on Animals and the National Institutes of Health (Protocol number: IS00000041-6).

No statistical methods were used to predetermine sample sizes for barcoding experiments, but sample sizes were similar to those reported in previous publications (LaFleur et al., 2019a, 2019b). Statistical methods were used to estimate sample sizes for screen experiments. Mice were randomized prior to pooling to obtain groups with comparable tumor sizes (mean and SD). Data exclusion was not used except for the barcoding and screen experiments, where tumors that did not grow beyond 50 mm³ in 10 days were excluded. Data collection and analysis were not performed blind to the conditions of the experiments.

The gRNA expression plasmid (pRDA_526) is available on Addgene (ID: 184859; Addgene). pRDA_526 is a modified form of our gRNA expression vector pXPR_071 (ID: 164558; Addgene) with the following changes: (1) modification of the tracrRNA, (2) modification of the stuffer between BsmBI digest sites, (3) addition of a gRNA capture sequence following the tracrRNA, and (4) replacement of the fluorophore Vex with Thy1.1. The pLX_305 vector (Broad Institute Genetic Perturbation Platform) was modified to express full-length OVA and is available on Addgene (ID: 184924; Addgene).

The gRNA oligonucleotides used in the 6,000-gRNA barcoding experiments were generated by sampling candidate gRNA sequences at NGG PAM (protospacer adjament motif) sites from the entire GRCm38 genome (October 2020 version) and eliminating sites that were inside any gene region (including introns/UTR or any site within 300 bases upstream of the transcription start site, according to the latest NCBI gene/transcript annotation at the time). In addition, gRNAs were eliminated if they had (1) any other matches to the genome where there was a 100% sequence match over the final 17 bases of the target 20mer plus the GG of the PAM, (2) an on-target score <0.5, and (3) any restriction sites or TTTT. This process was repeated until reaching 5,990. Ten Eef2-targeting gRNAs were included to bring the total library size to 6,000. Tables S1 and S3 contain the full list of gRNAs used in the 6,000-gRNA barcoding experiments.

The gRNA oligonucleotides used in the 100-gRNA screens were designed using the Broad CRISPR algorithm: https://portals.broadinstitute.org/gppx/crispick/public. We designed three gRNAs per gene for 30 genes and included five non-targeting control gRNAs and five intergenic-targeting control gRNAs. Table S4 contains the full list of gRNAs used in the 100-gRNA screens.

The UMI oligonucleotides used in the 6,000-gRNA barcoding experiments and the 100-gRNA screen experiments were generated by designing all permutations of oligonucleotides 6 bp in length with the restriction that each 6-bp UMI be at least 2 bp away from any other UMI (hamming distance = 2). Tables S1 and S3 contain the full list of UMIs used in the 6,000-gRNA barcoding experiments.

pRDA_526 was digested with BsmBI (Cat# R0580L; NEB). gRNA libraries and UMIs were simultaneously cloned into digested pRDA_526 via a three-part pooled ligation. Ligation products were cleaned up with isopropanol precipitation. Ligation products were then electroporated into Stbl4 bacteria (Cat# 11635-018; Thermo Fisher Scientific) and plated on large bioassay plates at 37°C. Plasmid libraries were then extracted using a HiSpeed Plasmid Maxi kit (Cat# 12662; Qiagen), quantified via NanoDrop, and sequenced to confirm even representation of the gRNAs. Plasmid libraries used were as follows: CP1656 (6,000-gRNA library with 96 UMIs/gRNA), CP1923 (6,000-gRNA library with 384 UMIs/gRNA), and CP1696 (100-gRNA library with 96 UMIs/gRNA).

293x cells (HEK variant obtained from Dr. Cigall Kadoch, Dana Farber Cancer Institute, Boston, MA, USA) were cultured in DMEM (Cat# 11995073; Gibco) supplemented with 10% FBS (Cat# F2442; Sigma-Aldrich) and 1% penicillin/streptomycin (Cat# 97063-708; VWR). B16-OVA cells were produced by transducing B16.F10 cells (a gift from Dr. Glenn Dranoff, Dana Farber Cancer Institute, Boston, MA, USA) with the pLX_305_Ovalbumin vector (Miller et al., 2019). B16-OVA cells were selected for 2 days in 2 µg/ml puromycin (Cat# 10766-886; VWR). 293x cells were not independently authenticated and B16-OVA cells were authenticated by flow cytometric analysis of SIINFEKL expression. Both cell lines were confirmed to be mycoplasma-negative and free of standard pathogens.

293x cells were transfected using LT-1 (Cat# MIR2305; Mirus Bio) with the packaging plasmids pCag-Eco (ID: 35617; Addgene), pMD2.G (ID: 12259; Addgene), and psPAX2 (ID: 12260; Addgene), and pRDA_526 plasmid containing the gRNA libraries. pCag-Eco and pMD2.G were used at a 1:1 mass ratio to make pseudotyped lentivirus. Lentivirus-containing supernatant was collected 48 h after transfection and spun at 480g for 5 min to remove cells. The supernatant was then frozen at −80°C until use.

Non-tissue culture–treated 6-well plates were coated with RetroNectin (Cat# T100B; Takara Bio) at ∼15 µg/cm² at 4°C overnight. The next morning, the RetroNectin was removed and plates were blocked with 2 ml/well of 2% BSA (Cat# A2153-50G; Sigma-Aldrich) in PBS^−/− (Cat# 82020-066; VWR) for 30 min at room temperature (22°C). The block was then removed and 2 ml of lentivirus was plated on RetroNectin-coated plates and centrifuged at 2,000g for 2 h at 32°C. Naive CD8⁺ T cells were isolated using a Naive CD8a⁺ T cell Isolation Kit (Cat# 130-096-543; Miltenyi) and then plated on virus-coated plates with 1 mg/ml LentiBOOST (Sirion Biotech) in complete RPMI media (RPMI 1640 with L-glutamine and no sodium pyruvate, pH 7.0–7.4 [Cat# 11-875-093; Thermo Fisher Scientific], 10% FBS, 1% penicillin/streptomycin, 10 mM HEPES [Cat# 97064-360; VWR], 55 µM 2-mercaptoethanol [Cat# 97064-588; VWR], 1% sodium pyruvate [Cat# 11360070; Thermo Fisher Scientific], and 1% non-essential amino acids [Cat# 11-140-050; Thermo Fisher Scientific]) containing 1 ng/ml recombinant murine IL-7 (Cat# 10780-184; VWR). Cells were cultured for 7 days in complete RPMI media with 1 ng/ml recombinant murine IL-7 to enable expression of Thy1.1 while keeping the cells naive (LaFleur et al., 2023). Media was changed on days 2 and 5 after transduction and fresh cytokines were added on day 4 after transduction. Flow cytometry was performed on day 7 after transduction to ensure the virus led to transduction, measured by Thy1.1 expression.

Samples were stained for 30 min on ice in the dark with antibodies in FACS buffer, containing PBS^−/− with 2 mM EDTA (Cat# 45001-122; VWR) and 1% FBS. Following staining, cells were washed twice with FACS buffer. Cells were then resuspended in FACS buffer and either (1) analyzed on a BD Celesta or (2) sorted on a BD Aria IIu. FACSDIVA Versions 9.0.1 and 9.4 were used for data collection. Compensation was performed with single-color controls. Gating was performed using fluorescence-minus-one controls where necessary. Samples were analyzed using Flowjo Version 10.8.1 software.

Biolegend: B220 (Clone RA3-6B2, Cat# 103208), CD25 (Clone PC61, Cat# 102012), CD44 (Clone IM7, Cat# 103006, 103031), CD45.1 (Clone A20, Cat# 110708), CD45.2 (Clone 104, Cat# 109824), CD62L (Clone MEL-14, Cat# 104417), CD8β (Clone YTS156.7.7, Cat# 126610, 126613, 126605), Gr-1 (Clone RB6-8C5, Cat# 108408), PD-1 (Clone 29F.1A12, Cat# 135214), Ter-119 (Clone TER-119, Cat# 116208), Thy1.1 (Clone OX-7, Cat# 202529), TruStain fcX (Clone 93, Cat# 101320), TCR Vα2 (Clone B20.1, Cat# 127822), and 7-AAD (Cat# 420404). Thermo Fisher Scientific: Near-IR Fixable Live/Dead (Cat# L34976).

Mice were anesthetized with 2.5% 2,2,2-Tribromoethanol (also known as Avertin) (Cat# AAA18706-22; VWR), shaved, and injected midway between the stomach and the flank subcutaneously with 5 × 10⁶ B16-OVA tumor cells. Tumors were measured every 2–3 days once palpable using a caliper. Tumor volume was determined by the volume formula for an ellipsoid: 1/2 × D × d², where D is the longer diameter and d is the shorter diameter. Mice were sacrificed when tumors reached 2 cm³ or upon ulceration, as per Institutional Animal Care and Use Committee guidelines.

Tumors were excised and mechanically minced. Tumors were then incubated in 2 mg/ml Collagenase Type 1 (Cat# LS004194; Worthington Biochemical) for 20 min at 37°C with gentle agitation. Collagenase was then neutralized with FBS-containing RPMI media (Cat# R8758-500ML; Sigma-Aldrich). Tumors were mechanically dissociated through a 70-µm filter. Lymphocytes were enriched using an Optiprep gradient (Cat# D1556; Sigma-Aldrich). CD8⁺ T cells were further enriched using CD8α microbeads (Cat# 130-117-044; Miltenyi).

LNs were isolated and mechanically dissociated through a 70-µm filter. For the in vivo screens, CD8⁺ T cells were then enriched using CD8α microbeads.

The 576,000-gRNA-UMI barcoding experiment was performed by transducing OT-1 CD45.1⁺/CD45.2⁺ Cas9-expressing naive CD8⁺ T cells with lentiviral supernatants from the CP1656 gRNA library, which contains a Thy1.1 reporter to identify transduced T cells. Naive CD8⁺ T cells were cultured for 7 days in recombinant murine IL-7. Naive CD8⁺ T cells were then stained with 7-AAD, Lineage (Ter119/Gr1/B220)-PE, and Thy1.1-BV421 and sorted (to >90% purity) for 7-AAD^– Lineage^– Thy1.1⁺ cells. Three sets of 500,000 cells were frozen at −80°C as input samples. We also sorted 3 × 10⁶ additional cells and transferred 300,000 cells intravenously to each of the 10 Cas9-expressing mice. The next day, these mice were injected subcutaneously with 5 × 10⁶ B16-OVA cells. Tumor growth was monitored and, when tumors reached ∼400 mm³, mice were euthanized and tdLNs and tumors individually processed as above. CD8⁺ T cell–enriched fractions from tumors and non-enriched LN samples were then stained and sorted (to >90% purity) for Near-IR Fixable Live/Dead^– CD8β⁺ CD45.1⁺/.2⁺ Thy1.1⁺ CD44⁺ TCR Vα2⁺ cells. Cell pellets were frozen down at −80°C. Genomic DNA was isolated from the three input samples, 10 tdLN samples, and 10 tumor samples using a DNeasy Blood and Tissue Kit (Cat# 69504; Qiagen).

The 2,304,000-gRNA-UMI barcoding experiment was performed by transducing OT-1 Cas9-expressing naive CD8⁺ T cells with lentiviral supernatants from the CP1923 gRNA library, which contains a Thy1.1 reporter to identify transduced T cells. Naive CD8⁺ T cells were cultured for 7 days in recombinant murine IL-7. To isolate Thy1.1⁺ T cells, naive CD8⁺ T cells were stained with Thy1.1-PE and then enriched on LS columns (Cat# 130-042-401; Miltenyi) using Ultrapure anti-PE-microbeads (Cat# 130-105-639; Miltenyi). 2.5 million Thy1.1⁺ CD8⁺ T cells were frozen at −80°C as an input sample. We then transferred a range of doses (6,000, 30,000, 60,000, 150,000, and 300,000) of Thy1.1⁺ CD8⁺ T cells intravenously to Cas9-expressing mice. The next day, these mice were injected subcutaneously with 5 × 10⁶ B16-OVA cells. Tumor growth was monitored and sets of mice were euthanized 6, 10, and 14 days after tumor implantation. CD8⁺ T cell–enriched fractions from tumors and non-enriched LN samples were then stained and sorted (to >90% purity) for Near-IR Fixable Live/Dead^– CD8β⁺ Thy1.1⁺ CD44⁺ TCR Vα2⁺ cells (LN) or Near-IR Fixable Live/Dead^– CD8β⁺ Thy1.1⁺ CD44⁺ cells (tumor). Cell pellets were frozen down at −80°C. Genomic DNA was isolated from the one input sample, 73 tdLN samples, and 73 tumor samples using a DNeasy Blood and Tissue Kit.

The in vivo CRISPR screen examining responses in the tdLN and tumor was performed by transducing a mix of CD45.1⁺/CD45.2⁺ and CD45.1⁺ OT-1 Cas9-expressing naive CD8⁺ T cells with lentiviral supernatants from the CP1696 library, which contains a Thy1.1 reporter to identify transduced T cells. Naive CD8⁺ T cells were cultured for 7 days in recombinant murine IL-7 (Cat #10780-184; VWR). Naive CD8⁺ T cells were then stained with 7-AAD, Lineage (Ter119/Gr1/B220)-PE, and Thy1.1-BV421 and sorted for 7-AAD^– Lineage^– Thy1.1⁺ cells. Four sets of cells containing 200,000, 643,000, 900,000, or 501,000 Thy1.1⁺ CD8⁺ T cells each were frozen at −80°C as input samples. We also sorted 16 × 10⁶ additional cells and transferred 300,000 intravenously to each of the 53 Cas9-expressing mice. The next day, these mice were injected subcutaneously with 5 × 10⁶ B16-OVA cells. Tumor growth was monitored and, when tumors began to regress, mice were euthanized. Tumors from 5 of the 53 mice did not grow to sufficient size or ulcerated and were excluded from analysis. On day 11, mice were pooled into six groups of eight mice per group for processing (grouping all eight tdLN or all eight tumors). Pooling was done by grouping the mice such that the pools had approximately equivalent average tumor size and SD. TdLNs and tumors were processed as above. CD8⁺ T cell–enriched fractions were then stained and sorted (to >90% purity) for Near-IR Fixable Live/Dead^– CD8β⁺ CD45.1⁺ (containing both CD45.1⁺ and CD45.1⁺/.2⁺) Thy1.1⁺ CD44⁺ TCR Vα2⁺ cells. Cell pellets were frozen down at −80°C. Genomic DNA was isolated from the four input samples, six tdLN samples, and six tumor samples using a DNeasy Blood and Tissue Kit.

The in vivo CRISPR screen examining multiple LNs was performed by transducing a mix of CD45.1⁺/CD45.2⁺ and CD45.2⁺ OT-1 Cas9-expressing naive CD8⁺ T cells with lentiviral supernatants from the CP1696 library. Naive CD8⁺ T cells were cultured for 7 days in recombinant murine IL-7 (Cat #10780-184; VWR). To isolate Thy1.1⁺ T cells, naive CD8⁺ T cells were stained with Thy1.1-PE and then enriched on LS columns (Cat# 130-042-401; Miltenyi) using Ultrapure anti-PE-microbeads (Cat# 130-105-639; Miltenyi). Two sets (both 100,000) of Thy1.1⁺ CD8⁺ T cells were frozen at −80°C as input samples. We then transferred 300,000 Thy1.1⁺ CD8⁺ T cells intravenously to each of the nine Cas9-expressing mice. The next day, these mice were injected subcutaneously with 5 × 10⁶ B16-OVA cells. Tumor growth was monitored and, when tumors began to regress, mice were euthanized. From each mouse, the ipsilateral inguinal LN, the ipsilateral axillary LN, the contralateral inguinal LN, the contralateral axillary LN, and the tumor were isolated. During processing, the nine mice were pooled into one group of nine mice for processing (grouping all nine ipsilateral inguinal LNs, all nine ipsilateral axillary LNs, all nine contralateral inguinal LNs, all nine contralateral axillary LNs, or all nine tumors). The LNs and tumors were processed as above. CD8⁺ T cell–enriched fractions were then stained and sorted for Near-IR Fixable Live/Dead^– CD8β⁺ Thy1.1⁺ CD44⁺ cells. Cell pellets were frozen at −80°C. Genomic DNA was isolated from the two input samples, four LN samples and one tumor sample using a DNeasy Blood and Tissue Kit.

We assessed DNA concentration using the Qubit dsDNA HS assay kit (Cat# Q32851; Thermo Fisher Scientific). We then performed a PCR using Titanium Taq (Cat# 638517; Takara) on 2 µl of some of the samples to ensure the DNA was of sufficient quality for amplification. Amplified samples were run on a 1.5% agarose gel to assess PCR product formation. The following primers were used for the PCR:

• 5′-AAT​GAT​ACG​GCG​ACC​ACC​GAG​ATC​TACA​CTC​TTT​CCC​TAC​ACG​ACG​CTC​TTC​CGA​TCTtcg​att​tct​tgg​ctt​tat​ata​tct​tgt​g-3′

• 5′-CAA​GCA​GAA​GAC​GGC​ATA​CGA​GATGTG​ACT​GGA​GTT​CAG​ACG​TGT​GCT​CTT​CCG​ATC​Tacc​gac​tcg​gtg​cca​ctt​ttt​caa-3′

Uppercase indicates P5/P7 flow cell attachment sequence; underline indicates Illumina sequencing primer; lowercase indicates vector primer binding sequence.

The Broad Institute Genetic Perturbation Platform amplified and indexed genomic DNA by PCR using Titanium Taq. A portion of the samples was then run on a gel to estimate product amounts among the samples. The samples were then pooled based on these estimates to achieve roughly equimolar representation in the pool. Samples were then sequenced on a HiSeq2000 with a 50-cycle kit (50-bp read 1, 8-bp index read).

Fastq files were demultiplexed using the Broad Institute Genetic Perturbation Platform’s program PoolQ (Version 3.3.2) into individual samples (sample indices are in Tables S1 and S3) with read counts per gRNA or gRNA-UMI pair. PoolQ source code is located on Github: https://github.com/broadinstitute/poolq. The 50-bp read is structured as follows:

• 5′-TGT​GTA​TCA​TAT​ATC​TCA​CGGTT​TGA​GAG​CTACCG​TTATAG​CAA​GTT​CAA-3′

Uppercase indicates gRNA sequence, underline indicates conserved pRDA_526 vector sequence, and bold indicates UMI sequence.

Upon adoptive transfer of cells, we assumed that two independent processes occur: first, many T cell clones underwent spontaneous extinction (“bottlenecking”). Second, many T cell clones expanded, in response to the antigenic challenge. We defined “clonal engraftment” as the fraction of transferred T cells (assumed to be clonally distinct) that survive initial bottlenecking. We sought to measure the clonal engraftment of transferred T cells based on an output sample (e.g., tumor or LN). The distribution of cells belonging to particular clones in this output sample was assessed via PCR of genomic DNA from the output sample using primers specific to the cells’ clone-identifying barcode, as described above. While the number of sequencing reads corresponding to a particular post-PCR barcode should be proportional to the number of input cells with that barcode, that number may be skewed by the stochastic nature of PCR (“PCR jackpotting”). We were interested in measuring the extent of bottlenecking, ignoring the effects of clonal expansion. To address this, we wanted to measure the probable number of T cell clones in the output sample, adjusting for the fact that multiple T cell clones may share a particular barcode.

To accomplish this, we assumed that each possible gRNA-UMI barcode has an equal likelihood of having a corresponding cell. For each possible gRNA-UMI barcode, the distribution of corresponding cell numbers k will follow a binomial distribution B(k, n, p), where n is the number of cells, and p = 1/(number of barcodes). In the limit of large n and small p, this may be approximated by the Poisson distribution Pois(k,λ = n*p).

We then estimated clonal engraftment as $Ncnumber of input clones$⁠. In the 6,000-gRNA, 96-UMI barcoding experiment, the number of input clones was assumed to be 300,000 and the number of possible barcodes (N_pb) was 576,000.

In developing FAUST, we examined multiple other existing algorithms developed to analyze either RNAi or CRISPR screening assays, namely RIGER (Luo et al., 2008), MAGeCK (Li et al., 2014), MAGeCKiBAR (Zhu et al., 2019), ZFC (Xu et al., 2021), BAGEL (Hart and Moffat, 2016), casTLE (Morgens et al., 2016), and STARS (Doench et al., 2016).

Broadly speaking, these techniques differ based on:

• ⋅How outliers—due to variability in construct efficacy, incomplete (e.g., heterozygous) gene KO, or PCR jackpotting—are accounted for

• ⋅How UMIs are utilized

• ⋅Whether the algorithm identifies gene essentiality or enrichment

• ⋅How the algorithm uses non-targeting constructs as controls

We found that all seven of the above algorithms ultimately calculate significance based on the likelihood of the null hypothesis that a given genetic target was associated with a test statistic that was distributed identically to the test statistics associated with all other targets in the screen. Three of the algorithms (MAGeCK, MAGeCKiBAR, and ZFC) all use different test statistics, but ultimately convert these to ranks (after ranking genetic targets according to the value of their corresponding summary statistic) and compute significance for each target using a modified version of Robust Rank Aggregation (Kolde et al., 2012). This approach (comparing each target to all other targets) is contrasted with FAUST, where each target is compared with an explicit set of controls. Specifically, FAUST tests the null hypothesis that the enrichment of gRNA-UMIs associated with a particular target is sampled from the same distribution as those of gRNA-UMIs associated with non-targeting or intergenic-targeting controls. The resultant significance of a target’s enrichment or depletion is therefore independent of the other targets included in a panel, except to the extent that such a “bystander” target may “out-compete” both targets and controls and reduce the statistical power to test the aforementioned null hypothesis. We found two algorithms (MAGeCKiBAR, ZFC) to be most applicable to FITS, and we, therefore, compared results from FAUST to those from these two algorithms. These algorithms were run with the wrapper function faust.utilities.run_alternative_test, with argument “test” set to be “mageckibar” or “zfc.”

Fold changes between the conditions were calculated by subtracting the log₂-normalized input gRNA counts from the tdLN or tumor log₂-normalized gRNA counts (for tdLN versus input and tumor versus input, respectively) or by subtracting the log₂-normalized tdLN gRNA counts from the log₂-normalized tumor gRNA counts in matched mice.

Fold changes between the conditions were calculated by subtracting the log₂-normalized input gRNA-UMI counts from the tumor log₂-normalized gRNA-UMI counts.

Heatmaps were generated using the seaborn v0.12.0 seaborn.heatmap subroutine. Cluster maps (hierarchically clustered heatmaps) were generated using the seaborn v0.12.0 seaborn.clustermap subroutine.

Briefly, log₂-normalized gRNA count tables for input, replicate pools of mice from the tumor, and replicate pools of mice from the tdLN were created as above. TdLN vs. input, tumor vs. input, and tdLN vs. tumor fold changes were calculated, as above, for replicate mouse pools A-F. The Pearson correlation coefficient was calculated for all pairwise combinations of A-F for a given comparison (tdLN vs. input, tumor vs. input, or tdLN vs. tumor). Each pairwise combination is indicated by an individual square in the heatmap.

First, a designated number of UMIs were chosen to be excluded (all values, inclusive, between 0 and 95). gRNA-UMI counts corresponding to these UMIs were removed. Then, counts for all gRNAs were summed together. Counts were then normalized on a per-sample basis such that the sum of gRNA counts for all gRNAs for a given sample would be 1. We then computed the fold change between these normalized counts between tumor and input (using the normalized counts of the mean of three input aliquots), LN and input (input as before), and tumor and LN (where tumor and lymph were sample-matched). We then computed the ICC2 using the pingouin (https://github.com/raphaelvallat/pingouin, v0.5.3) function, pingouin.intraclass_corr. Therein, “targets” represented the different gRNA, “raters” the different biological replicates, and “ratings” the values of the fold change. This procedure was repeated 1,000 times, with randomly selected (using np.random.choice) UMIs upon each iteration. We took the average intraclass correlation over these 1,000 iterations. Plots show mean values for ICC2, i.e., “Single random raters.”

FAUST is an open-source Python algorithm available at https://github.com/scyrusm/faust.

Outputs from PoolQ were aggregated with the function faust.utilities.read_gpp_output, which combines UMI-specific files into a single dataframe, with each line corresponding to the site-specific counts for each gRNA-UMI. The function faust.utilities.get_summary_df was then used to compute P values and effect sizes. Sets of gRNA-UMIs were designated as “controls” (corresponding to gRNA that target non-coding regions of the genome or that target no regions of the genome) and “experimental” (corresponding to gRNA targeting a specific gene), and sites were designated as “input” and “output.” In the case of overall proliferation, “input” represented the preinjection counts of gRNA-UMI and “output” those in a tumor or LN. To study multiple sites, “input” represented the source site (e.g., an LN) and “output” the destination (e.g., a tumor). faust.utilities.get_summary_df was then used to compute ratios R_i,jk for each gRNA-UMI i, for input j and output k. Where a count threshold was used, any count value less than or equal to the threshold was regarded as zero. Any undefined ratios (for example, if there were zero detected counts for a particular gRNA-UMI in the designated input) were not included in downstream calculations. The Mann–Whitney U statistic and associated P values were computed (via the scipy function scipy.stats.mannwhitneyu, by default using optional argument “alternative” set to “two-sided”) between the set of R_i,jk for i in S_c, i’ in S_e, where S_c and S_e represent the set of control and experimental gRNA-UMIs, respectively. The common language effect size was then computed as U/(|S_c|*|S_e|), where |S| denotes the number of considered (i.e., finite) elements in S.

The Benjamini–Hochberg q-value corresponding to each P value was then computed using the routine faust.utilities.nan_fdrcorrection_q.

Where there were multiple replicates in a single screen, aggregation thereof was performed in one of two modes: liberal and conservative. In the liberal setting, we tested the intersection null hypothesis i.e., the hypothesis that “the null hypothesis (that the median of ratios between input and output is the same for control guides as for experimental guides with a particular gene target) holds in all replicates.” This intersection null is evaluated via Fisher’s combined test, with P values for each replicate calculated as above (via Mann–Whitney U). In the conservative setting, we tested the union null hypothesis, i.e., the hypothesis that “the null hypothesis (that the median of ratios between input and output is the same for control guides as for experimental guides with a particular gene target) holds in at least one replicate.” To test the union null hypothesis, we used the maximum P value over all replicates.

To determine if the 6,000-gRNA library and the 96 UMIs present in the barcoding experiment impacted cell fitness, we analyzed the tdLNs and tumors across the 10 replicate animals using FAUST. To assess the impact of particular gRNAs, we compared either the ratio (input versus tdLN, or input versus tumor) distribution of the 96 gRNA-UMIs corresponding to a given gRNA versus all other gRNA-UMIs (downsampled randomly to 0.5% of gRNA-UMIs, i.e., 2,880 total) and computed effect sizes for each gRNA per the FAUST algorithm. To assess the impact of particular UMIs, we compared the ratio (input versus tdLN, or input versus tumor) distribution of the 6,000 gRNA-UMIs corresponding to a given UMI versus the 570,000 other gRNA-UMIs with different UMIs. Results were summarized across 10 replicates by taking the median across replicates of the common language effect size of the comparisons described above.

The number of cells with a particular gRNA (or genetic target) was calculated under a Poisson approximation, which is assumed to be valid when the number of cells is much smaller than the number of UMIs. This was estimated as N_total*(ln(N_total) − ln(N_total − N_detected)), where N_detected is the number of UMIs corresponding to a particular gRNA or target with counts above the threshold and N_total is the total number of UMIs. This represents the sum $∑x=1∞$ Poisson(λ)(x), for λ = N_detected/N_total. This was implemented in the routine faust.utilities.estimate_cell_input.

Swarm/violin plots displaying top hits from the FAUST algorithm were generated using the routine faust.visualization.plot_top_hits.

To assess how well in vivo screening results were modeled by a negative binomial distribution, we used the R (v4.0.5) package “vcd” (v1.4.10) and the contained function “goodfit” to calculate optimal negative binomial distribution according to either a maximum likelihood or a minimum χ² statistic. For each of these criteria, we compared the theoretical (based on the fit negative binomial hyperparameters) and empirical moments (mean, SDs, skews, and kurtosis) as calculated using the functions numpy.mean, np.std (numpy v1.21.6), scipy.stats.skew, and scipy.stats.kurtosis (scipy v1.9.3). To test whether the median of theoretical values was different from that of empirical, we applied a Wilcoxon signed-rank test using the function scipy.stats.wilcoxon.

Cell coverage was estimated by taking the engraftment rate of 2% and a cell transfer dosage of 300,000 cells, which equates to ∼6,000 clones engrafting each mouse. In the 100-gRNA screen (Fig. 2), 8 mice were pooled, equating to ∼48,000 clones after pooling. The library has 9,600 gRNA-UMIs, so we estimate each gRNA-UMI is represented by ∼5 cells.

We noted that barcodes with high representation for a given UMI had higher-than-expected representation for other UMIs, violating the assumption that each gRNA-UMI is clonally distinct and uncorrelated (particularly so for day 14 tumor samples). This is likely due to read errors in the UMI portion of the read (UMIs differ by as few as two base pairs). To address this, we implemented the following models:

M1: let n_recovered be the number of recovered cells from a sample and N be the total number of distinct gRNA-UMI (here 384*6000 = 2.304*10⁶). Assuming that each of the recovered cells is clonally distinct, the number of gRNA-UMIs with non-zero count will be approximately Poisson distributed with rate λ_true = n_recovered/N. Then the likelihood of observing j distinct gRNA-UMIs, conditional on having observed one or more, is $Ptrue(j,λtrue)=1−Pois(λtrue)(j)1−Pois(λtrue)(1)$⁠, where Pois(λ_true) is the cumulative distribution function of the Poisson distribution with rate λ_true. We additionally assumed a fixed read error rate of 1e-3, amounting to an approximately 1e-6 two-read error rate, which is the minimum distance between any two UMIs in the study. Assuming this, the likelihood of there being a gRNA-UMI that is in fact a read error of the gRNA-UMI with the same barcode but different UMI is approximately P_error(λ_error) = 1 − Pois(λ_error)(1), where λ_error = n_highest* 1e-6, where n_highest is the number of reads assigned to the gRNA-UMI with the highest number of reads within a set of gRNA-UMIs with a particular barcode (here we assume that the reads stemming from read errors are Poisson distributed with rate λ_error across the 384 UMIs). To select a read threshold where the likelihood of observing j distinct gRNA-UMIs exceeds that of observing one read error, we identified the lowest integer j such that P_true(j, λ_true) < P_error(λ_error). We then ranked the counts of reads for each gRNA-UMI corresponding to a particular barcode and let the j^th count represent a barcode-specific read threshold, whereby all counts with that many or fewer reads were reassigned zero reads.

We furthermore noticed that some samples had higher apparent engraftment than the number of cells recovered (particularly for day 6 tumors). This appeared to be due to sample-level mis-assignment, possibly due to read errors in the index sequence. To address this, we implemented the following model:

M2: For a given sample, for each gRNA-UMI, let counts <1e-5 * n_{max_nonsubject} be set to zero, where n_{max_nonsubject} is the maximum number of counts for that gRNA-UMI for all other samples not from that subject (i.e., not corresponding to two tissues, e.g., LN and tumor, from a single biological subject, or to the input aliquot).

We jointly incorporated these models by regarding all gRNA-UMI counts as zero where they do not pass the filtering requirements of either M1 or M2.

The “expansion factor” (EF), as used in Fig. S5 i, was computed as follows: EF at input was set to 1 for all tissues. EF at LN (EF_LN) was set to be E_CLES/(1 − E_CLES), where E_CLES is the common language effect size of a gene in question compared to control, as computed by FAUST using input-to-LN comparison. EF at tumor via a particular LN (⁠$EFLN→T$⁠) was set to be EF_LN (E_CLES/(1 − E_CLES)), where here E_CLES is the effect size of a gene in question relative to control, as computed by FAUST using LN-to-tumor comparison.

All code and scripts used for data analysis related to the FAUST algorithm are available at https://github.com/scyrusm/faust. Other code notebooks used in the development of this manuscript are available at https://github.com/scyrusm/fits.

Fig. S1 shows characterization data of transduced naive CD8⁺ T cells and analyses of barcoding read count thresholding. Fig. S2 shows analyses of gRNAs and UMIs in the 6,000-gRNA barcoding experiment with 96 UMIs/gRNA as well as tumor curves and cell recovery from the 6,000-gRNA barcoding experiment with 384 UMIs/gRNA. Fig. S3 shows tumor growth curves, cell recovery, gRNA-UMIs detected, gRNAs detected, and tdLN versus input fold change plots from the 100-gRNA screen. Fig. S4 shows the distribution of sequencing reads relative to the binomial distribution, hit calling using the liberal mode of FAUST from the 100-gRNA screen, and comparisons of FAUST to MAGeCK for the 100-gRNA screen. Fig. S5 shows tumor growth curves, cell recovery, and gRNA-UMI recovery from the multisite 100-gRNA screen as well as analyses of gRNA-UMI clonality. Table S1 includes the gRNAs, UMIs, and sample indices used in the 6,000-gRNA barcoding experiment with 96 UMIs/gRNA. Table S2 includes the gRNA and gRNA-UMI count data from the 6,000-gRNA barcoding experiment with 96 UMIs/gRNA. Table S3 includes the gRNAs, UMIs, and sample indices used in the 6,000-gRNA barcoding experiment with 384 UMIs/gRNA. Table S4 includes the gRNAs, UMIs, and sample indices used in the 100-gRNA screens. Table S5 includes the gRNA counts, the gRNA-UMI counts, and hit calling for the 100-gRNA screen. Table S6 includes FDR data for MAGeCKiBAR and ZFC for the 100-gRNA screen. Table S7 includes the gRNA counts, the gRNA-UMI counts, hit calling, and MOI for the 100-gRNA screen experiment with multiple LNs.

All barcoding and screen sequencing data have been deposited on NCBI’s Gene Sequence Read Archive (BioProject: PRJNA957312). The other data and materials that support the findings of this study are available from the corresponding author upon reasonable request.

We thank Thomas Green, Xiaoping Yang, and David Root from the Broad Institute Genetic Perturbation Platform for assistance with barcoding and screen experiment design and analysis.

This work was supported by funding from grant no. U19AI133524 from the National Institute of Allergy and Infectious Diseases to A.H. Sharpe and J.G. Doench and sponsored research funding from Merck Sharp & Dohme LLC, a subsidiary of Merck & Co., Inc., Rahway, NJ, USA. S.C. Markson was supported through the American Association of Immunologists Intersect Fellowship Program for Computational Scientists and Immunologists.

Author contributions: M.W. LaFleur and A.H. Sharpe conceived the project. L.E. Milling, N.M. Derosia, G.H. Hickok, T.H. Nguyen, and M.W. LaFleur designed experiments. L.E. Milling, Q. Tjokrosurjo, N.M. Derosia, I.S.L. Streeter, G.H. Hickok, A.M. Lemmen, T.H. Nguyen, P. Prathima, and M.W. LaFleur acquired data. S.C. Markson, Q. Tjokrosurjo, and W. Fithian designed and performed analyses of the barcoding and screen experiments. M.A. Schwartz and N. Hacohen assisted with lentiviral transduction approaches. J.G. Doench assisted with screen methodology and analyses. L.E. Milling, S.C. Markson, Q. Tjokrosurjo, M.W. LaFleur, and A.H. Sharpe wrote the manuscript. All authors edited the manuscript. M.W. LaFleur and A.H. Sharpe supervised the project.