Mutational processes shape the landscape of TP53 mutations in human cancer.

Unlike most tumor suppressor genes, the most common genetic alterations in tumor protein p53 (TP53) are missense mutations1,2. Mutant p53 protein is often abundantly expressed in cancers and specific allelic variants exhibit dominant-negative or gain-of-function activities in experimental models3-8. To gain a systematic view of p53 function, we interrogated loss-of-function screens conducted in hundreds of human cancer cell lines and performed TP53 saturation mutagenesis screens in an isogenic pair of TP53 wild-type and null cell lines. We found that loss or dominant-negative inhibition of wild-type p53 function reliably enhanced cellular fitness. By integrating these data with the Catalog of Somatic Mutations in Cancer (COSMIC) mutational signatures database9,10, we developed a statistical model that describes the TP53 mutational spectrum as a function of the baseline probability of acquiring each mutation and the fitness advantage conferred by attenuation of p53 activity. Collectively, these observations show that widely-acting and tissue-specific mutational processes combine with phenotypic selection to dictate the frequencies of recurrent TP53 mutations.

TP53 is the most frequently mutated gene in human cancer[1]. Through interrogation of genome-scale RNAi[11] and CRISPR-Cas9[12] loss-of-function screens, we found that cells harboring wild-type (WT) TP53 acquired a selective fitness advantage when TP53 itself or several upstream activators (TP53BP1, CHEK2, ATM, USP28) or effector genes (CDKN1A, ZMAT3) were suppressed or deleted (Fig. 1a, c, d, Supplementary Fig. 1a, and Supplementary Table 1). We also found that suppression or deletion of TP53 was well-tolerated in cells harboring TP53 missense mutations, indicating that gain-of-function (GOF) activities associated with these mutations[6-8] did not affect the propagation of cancer cell lines (Fig. 1a and Supplementary Fig. 1a). These observations suggest that loss of canonical p53 function plays a central role in the selective advantage associated with TP53 mutation. However, because these assays rely on natural TP53 variation present in cancer cell lines, they do not provide functional information about every possible TP53 mutation.

Figure 1

Deletion of endogenous wild-type but not mutant TP53 impacts fitness in human cancer cells

Comparison of enrichment scores for CRISPR-Cas9 reagents targeting TP53 (a) or MDM2 (b) in cell lines whose TP53 statuses were defined using TP53 mutation, copy-number, target gene expression, and nutlin-3 sensitivity data (see Supplementary Table 1). Only cell lines with concordant functional and genetic classifications were included in these analyses. Each point represents the gene-level score for a given cell line, and error bars indicate the mean and standard deviation of each group. Cell lines in the loss-of-function (LOF) category harbor homozygous deletion, frameshift, or nonsense mutation of TP53 and express low levels of p53 protein (see Supplementary Fig. 1c). (****, P < 0.0001, Two-tailed Welch’s t-test) (c) PARIS, a rescaled normalized mutual information (RNMI)-based statistical analysis[52] was used to nominate genes whose enrichment scores were significantly different between the p53 non-functional mutant (LOF, missense, splice-site, and indel, with functional score < 0) and p53 functional wild-type (functional score > 0) cell line classes in genome-scale CRISPR-Cas9[12] and RNAi[11] screens. Reported p53 pathway components that scored as significant (FDR < 0.05) in both analyses are highlighted in red. (d) Infographic depicting differential enrichment scores of reported p53 pathway members and target genes in genome-scale CRISPR-Cas9 (lower left triangle) and RNAi (upper right triangle) screens. Dashed lines indicate transcriptional regulation. Genes scoring as significant in both analyses are outlined in yellow.

To determine the function of each missense or nonsense TP53 mutation, we created a comprehensive library of TP53 mutants using Mutagenesis by Integrated TilEs (MITE)[13-15]. We designed the library such that each allele would contain a single mutation, and that each of the 20 natural amino acids and a stop codon would be represented at each codon position. We also included alleles harboring silent mutations, which did not change the amino acid sequence of p53 but allowed us to track WT allele performance. Our mutagenesis approach generated greater than 99.8% of the expected mutant alleles, as gauged by massively-parallel sequencing of the expression plasmid pool (Supplementary Table 2).

To study the function of these alleles in the presence or absence of endogenous p53, we created isogenic TP53-WT (p53WT) and -null (p53NULL) A549 human lung carcinoma cell populations using CRISPR-Cas9-mediated gene editing[16]. These cells display potent p53-dependent drug responses, proliferate robustly in culture, and allow for stable expression of lentivirally-delivered TP53 at near-endogenous levels when present at a single copy per cell (Supplementary Fig. 2a). We capitalized on the differential responses of these isogenic cells to two p53-activating agents, nutlin-3 and etoposide, and performed pooled positive-selection screens designed to enrich for dominant-negative (DN), loss-of-function (LOF), or WT-like alleles.

Nutlin-3 interrupts the interaction between p53 and MDM2[17], a key negative regulator of p53 stability[18] and transcriptional function[19]. Suppression of MDM2 activity is critical for oncogenic stress-[20-22] and DNA damage-induced[23] p53 activation, and cancer cell lines harboring WT p53 require MDM2 expression for proliferation (Fig. 1b–d and Supplementary Fig. 1b). As expected, nutlin-3 treatment impaired the proliferation of p53WT A549 cells but had no effect on p53NULL cells (Supplementary Fig. 2b–c). Expression of exogenous WT p53 restored nutlin-3 sensitivity in p53NULL cells but did not alter the nutlin-3 response in p53WT cells. Notably, expression of mutant p53 p.Pro278Ala had no effect on p53NULL cells but rendered p53WT cells partially nutlin-3-resistant, indicating that this allele is deficient for WT function, and also interferes with endogenously expressed WT p53 in a DN manner. In the heterozygous setting, similar transcriptionally-inactive full-length p53 mutants sequester WT p53 from active complexes, and thus reduce p53 function more efficiently than mutants that are unstable or prematurely truncated[24,25]. We anticipated that we could identify all such DN TP53 alleles by introducing the MITE library into p53WT cells and treating them with nutlin-3. Performing the complementary experiment in p53NULL A549 cells would allow for the selection of LOF alleles and the depletion of alleles with WT-like activity.

Using a competition assay, we found that p53NULL A549 cells were more sensitive than their p53WT counterparts to high doses of DNA double strand break (DSB)-inducing agents, most notably etoposide and doxorubicin (Supplementary Fig. 3). Although in mouse thymocytes, DNA damage-mediated induction of WT p53 leads to apoptosis[26-28], in other contexts, p53 activation induces cell cycle arrest and DNA damage repair and thereby antagonizes p53-independent modes of cell death that result from unresolved DNA damage[29-31]. Our observations indicate that in A549 cells, WT p53 promotes cell survival in response to high doses of DSB-inducing agents.

Indeed, we found that forced expression of WT p53 in p53NULL A549 cells prevented the dramatic cell death induced by 5 μM etoposide treatment, whereas expression of mutant p53 p.Pro278Ala or Renilla luciferase, a negative control, had no effect (Supplementary Fig. 2e). We anticipated that by introducing the p53 MITE library into p53NULL cells and treating them with etoposide, we could enrich for alleles bearing WT-like activity and deplete LOF alleles. We noted that although the p.Pro278Ala mutant exhibited LOF in p53NULL cells, it failed to interfere with the etoposide response in p53WT cells, suggesting that its DN activity is limited to a subset of p53 responses (Supplementary Fig. 2d).

To perform the enrichment screens, we introduced the p53 MITE library into p53WT and p53NULL cells at a low multiplicity of infection (MOI), selected cells expressing TP53 alleles using puromycin, and treated these populations with nutlin-3 or etoposide for 12 days (Fig. 2a). We then enumerated which TP53 alleles were present by massively-parallel sequencing. We defined allele enrichments using the log2-fold change in reads in experimental arms relative to early time point samples and identified alleles that were selected under each of these conditions (Fig. 2b–d). The resulting allele enrichment scores are shown as heat maps juxtaposed to the reported p53 domain structure and the frequency of somatic TP53 mutations in human tumors[2,32] (International Association for Research on Cancer (IARC) database R18; see URLs) (Fig. 2e). We found that alleles harboring silent mutations were significantly depleted in p53NULL cells treated with nutlin-3 and enriched in p53NULL cells treated with etoposide relative to alleles bearing premature stop codons, suggesting that these assays effectively discriminated between WT and null alleles (Fig. 2f–h, P < 0.0001, Wilcoxon rank sum). We also found that most alleles, especially those with mutations in the DNA-binding domain (DBD), performed similarly in all screen conditions (Supplementary Fig. 4). Indeed, of the 3819 DBD missense variants we tested, 1481 showed evidence of LOF in one or both assays performed in the absence of p53, and 1219 (82.3%) of these exhibited DN activity (Supplementary Fig. 5). In addition, we identified several mutations in the transactivation domain (positions 18, 19, 22, 23, 24, 25, 26, 27, 29, 30), DNA-binding domain (positions 276, 280, 281), and tetramerization domain (positions 330, 332, 337, 338, 341, 344, and 348) that exhibited differential responses in the three assays. We also found that not all alleles bearing stop codon mutations behaved as LOF alleles, with different patterns observed for stop codons at positions 1–43, 44–289, 290–362, and 363–393 (Fig. 2b–d, Supplementary Fig. 4–5). We also noted significant associations between the enrichment score of each allele, its transcriptional activity measured in yeast[33], and the degree of evolutionary conservation of the mutated residue, as assessed by Align-GVGD[34] (Spearman rank, P < 0.0001) (Supplementary Fig. 4 j–k). To facilitate interpretation of these datasets and to accommodate new datasets as they become available, we created a website that allows the function of any tested TP53 allele to be queried (PHenotypic Annotation of TP53 Mutations - PHANTM; see URLs). The raw read counts for each allele in each condition are provided in Supplementary Table 2 and mean Z-scores are reported in Supplementary Table 3. We note that the effect of missense mutations on splicing cannot be inferred from our data, and therefore mutations at residues adjacent to splice junctions (codons 25, 26, 32, 33, 125, 126, 187, 188, 224, 225, 261, 262, 307, 308, 331, 332, 367, and 368) may also alter TP53 through effects on splicing.

Figure 2

Comprehensive mutational scanning of TP53

(a) A library comprising 8258 mutant TP53 alleles was introduced into A549 p53WT and p53NULL cells in a pooled format under conditions that favor the integration of a single vector in each cell. Library-infected p53WT cells were treated with nutlin-3, and library-infected p53NULL cells were treated with either nutlin-3 or etoposide. After 12 days, genomic DNA was harvested, PCR-amplified and subjected to next generation sequencing. (b–d) Heat maps of normalized allele enrichment scores (Z-scores) with codon-level average Z-scores plotted at right. (e) Left, The reported domain structure of p53 with residues 175, 248, and 273 highlighted: TAD = transactivation domain, PRD = Proline-rich domain, DBD = DNA-binding domain, ZN = Zinc-binding domain, 4D = tetramerization domain, CTD = C-terminal domain. Right, Total number of missense and nonsense mutations found at each codon in the IARC database[2,32]. (f–h) Density plot of alleles with silent mutations (wild-type alleles), nonsense mutations at codons 44–289 (loss-of-function alleles), missense mutations that are common in cancer, and SNV-generated missense mutations that have never been observed in cancer. Differences among all groups of alleles were significant in each condition (P < 0.0001, Wilcoxon rank sum test).

The observation that 30% of all p53 missense mutations found in human tumors affect five residues (Arg273, Arg248, Arg175, Arg282, and Gly245)[2,35] suggests that these residues play key roles in p53 function or that these positions are selectively targeted by mutational processes in cancer. When we examined these hotspot missense alleles, we found that while they exhibited LOF and DN activity, they were not the top scoring alleles in the library (Supplementary Fig. 4l). Moreover, we failed to detect TP53 dependency in cancer cell lines harboring endogenous hotspot mutations (Supplementary Table 1). Thus, although these residues are important for p53 function, these observations argue that mutations may be common at these sites because of the inherent mutability of the DNA sequences that encode them. Indeed, all five of these codons contain methylated CpG dinucleotides, which are highly susceptible to mutation via spontaneous deamination of the cytosine residue[36], and correspondingly, most of the substitutions found at these positions are C>T mutations.

Recently, unbiased approaches to identify and quantify the activity of mutational processes in cancer have been developed using whole-genome and whole-exome data from more than 10,000 human tumor samples[9,10]. These efforts have uncovered 30 distinct mutational signatures that describe the likelihood of acquiring a specific base change, given the activity of an underlying mutational process and the identity of the bases that immediately flank the mutated base – a parameter termed the trinucleotide mutation context (Fig. 3a) (COSMIC Mutational Signatures; see URLs). Of the 30 mutational signatures that have been identified, 17 have been attributed specific etiologies[10], many of which have previously been associated with specific TP53 mutations[36-43] (Fig. 3b).

Figure 3

Tissue-of-origin-selective TP53 mutations are linked to specific mutational processes

(a) Each signature in the COSMIC Mutational Signatures database contains 96 mutation probabilities, one for each trinucleotide mutation type[9,10]. To assign a baseline mutation probability to each TP53 allele for each signature, we first determined its trinucleotide mutation type, and then assigned the corresponding value from the database. Depicted here is the assignment of a mutation probability for Arg248Gln under the influence of Signature 1. (b) Fisher’s exact tests were performed to identify TP53 mutations that occur significantly more frequently in specific tumor types in the IARC[2,32] and GENIE[35] databases. The heatmap shows the relative mutation probabilities for each of the indicated TP53 mutations under the influence of each signature, depicted graphically on a white-to-red scale as the signature-specific percentile of all SNV-derived TP53 mutations (n = 2810). The percent of tumors of each tissue-of-origin in which a given signature was found by Alexandrov et al.[10] is depicted graphically on a white-to-black scale.

We hypothesized that the tissue-selective activity of mutational processes would manifest as tissue-selective enrichment of certain TP53 mutations. To identify such mutations, we queried two independent TP53 mutation databases (IARC[2,32] and GENIE[35]) and found 25 mutations that were significantly over-represented in a specific tumor type in both databases (Fig 3b). We next assigned each TP53 mutation a baseline mutation probability as a function of each mutational signature, and found that 23 of the 25 tissue-selective TP53 mutations were associated with mutational signatures commonly found in the relevant tissue-of-origin. Among these were mutations enriched in colon, skin, bladder, lung, and liver cancers that are associated with signatures of mismatch repair deficiency (Signature 6), UV exposure (Signature 7), APOBEC activity (Signatures 2 and 13), tobacco smoke (Signature 4*), and aflatoxin exposure (Signature 24), respectively (Fig 3b). We also noted that the most frequent TP53 mutations found in all tumor types were associated with the most commonly-observed mutational signature, Signature 1, which reflects the natural decomposition of 5-methylcytosine to thymine associated with aging[9,10]. Additionally, we found that several TP53 alleles, which exhibited WT-like activity in our assays but are recurrently observed in human tumors, exhibit high intrinsic mutability, as judged by Signature 1, suggesting that this common mutational process likely generates recurrent TP53 passenger mutations (Supplementary Figs. 6–7). Collectively, these findings suggest that mutational processes play an important role in shaping the landscape of TP53 somatic driver and passenger mutations.

We therefore attempted to model the TP53 mutational spectrum as a function of (i) gene-agnostic mutational signatures[9,10] and (ii) phenotypic selection as assessed by our TP53 mutagenesis screen data. We included parameters for each of the three screen conditions and projections of mutational signatures that are found across all cancer types (Signatures 1 and 5) as well as those that drive tissue-selective TP53 mutation patterns (Signatures 2, 4*, 6, 7, 13, and 24) (Fig. 3b). The models were trained to predict the frequency of somatic TP53 mutations found in the IARC database[2,32] and were validated in the independent GENIE database[35]. We found that models derived using only mutational signatures or phenotype scores yielded weak but significant correlations with the number of mutations observed at each codon position in p53 in human tumors (Pearson R2 = 0.11 and 0.29 respectively, P < 0.0001). However, a combined model that included all parameters recapitulated the observed mutation frequency at each position with high accuracy, capturing both DBD and hotspot enrichment patterns (Pearson R2 = 0.77, P < 0.0001) (Fig. 4a). This combined model was not simply over-fitting the training dataset, as it performed equally well in the validation database (Pearson R2 = 0.77, P < 0.0001) (Fig. 4b). We also verified that the models showed similar performance when predicting the frequency of mutations at a test set of codon positions not used for model-fitting (Supplementary Fig. 6). Ultimately, this model suggests that mutational processes create genetic diversity at the TP53 locus in somatic tissues using the same rules that govern their activity genome-wide, and that cells acquiring LOF and DN TP53 mutations have a fitness advantage over those that retain WT p53 function. This model provides a parsimonious explanation for the enrichment of TP53 missense mutations generally and hotspot mutations, in particular, in human tumors.

Figure 4

The TP53 mutational spectrum modeled as a function of mutational signatures and phenotypic selection

(a) Generalized linear models were trained to predict the mutation frequency of each TP53 allele in the IARC database[2,32] using mutational signatures from the COSMIC database[9,10], phenotypic selection data from the TP53 MITE library screens, or both. (b) Position-level mutation rates predicted by the combined model are plotted downwards and observed mutation rates in the GENIE validation database are plotted upwards.

It has previously been argued that loss of p53 function is the critical determinant that underlies the selection of TP53 mutations in cancer[33]. However, the preponderance of missense mutations relative to truncation mutations also argues that full-length mutant p53 actively promotes tumor development. Through our comprehensive screening approach, we found that more than 80% of full-length p53 DBD missense mutants that display LOF also display DN activity, suggesting that the ability of mutant p53 to interfere with WT p53 is critical during tumorigenesis. Although other p53 GOF activities may affect different aspects of tumor biology not assessed in this study, the lack of mutant p53 dependency in cancer cell lines coupled with the robust relationships between LOF, DN activity, and mutation frequency indicate that the observed spectrum of TP53 mutations likely arises due to the selection of LOF and DN alleles that are generated by specific mutational processes.

Methods

TP53 annotation of human cancer cell lines

The functional and genetic TP53 status of 966 cell lines was determined using the Cancer Cell Line Encyclopedia (CCLE; see URLs)[44], Genomics of Drug Sensitivity in Cancer (GDSC; see URLs)[45], Cancer Target Discovery and Development (CTD2; see URLs)[46], and The Cancer Genome Atlas (TCGA accessed via cBioPortal; see URLs)[47] databases. Cell lines were first separated into two functional classes by considering nutlin-3 sensitivity data from GDSC and CTD2, and a p53 target gene expression signature[48] computed using CCLE data. Each cell line was provisionally considered as p53 functional if the functional score (calculated as [Target Genes CCLE Z-score] - [Nutlin-3 CTD2 Z-score] - [Nutlin-3 Sanger Z-score]) was above 0, and provisionally considered as p53 non-functional if this value was below 0. Cell lines in the p53 functional class were declared p53 wild-type (WT) if no TP53 alterations were detected by CCLE, GDSC, or TCGA (n = 252), and discarded as ambiguous if any TP53 alterations were found (n = 104). Cell lines in the p53 non-functional class were declared p53 mutant if any genetic TP53 alteration was found (n = 528) and discarded as ambiguous if no TP53 alterations were found (n = 82). The p53 mutant class was further divided into four subclasses: a loss-of-function (LOF) subclass, comprising cell lines with nonsense mutations, frameshift mutations, or homozygous deletions; a missense subclass; a splice-site subclass; and an in-frame insertion/deletion subclass. Cell lines with multiple TP53 alterations were classified using the following precedence order: missense > in-frame > splice site > LOF. Refer to Supplementary Table 1 for the full classification matrix.

Analyses of Project Achilles screening data

Pooled genome-scale LOF screens have been described[11,12]. The p53 classification scheme allowed us to annotate 348 cell lines in the RNAi dataset[11] (LOF, n = 67; missense, n = 150; splice-site, n = 21; in-frame, n = 7; WT, n = 103) and 276 cell lines in the CRISPR-Cas9 dataset (LOF, n = 45; missense, n = 130; splice-site, n = 14; in-frame, n = 4; WT, n = 83). Gene-level enrichment scores for RNAi screens were computed using DEMETER[11], an algorithm that maximizes the contributions of on-target reagents and controls for microRNA-like seed effects to reduce the contributions of off-target reagents. CERES[12], an algorithm that minimizes the contributions of copy-number effects was used to compute gene-level scores for the CRISPR-Cas9 screens. PARIS (Probability Analysis by Ranked Information Score), a permutation-based analytical tool was used to identify differential gene-level enrichment scores between cell lines in the p53 WT class and cell lines in the p53 mutant class (GenePattern; see URLs). A p53 Pathway Score was generated for each cell line using either DEMETER or CERES scores using the following formula: (MDM2 + MDM4 + PPM1D + USP7 + UBE2D3) – (TP53 + TP53BP1 + CHEK2 + ATM + CDKN1A).

Deletion of TP53 in A549 cells using CRISPR-Cas9

A549 cells were seeded in 6-well dishes (Costar, Corning, NY, USA) in normal culture media (DMEM containing 10% FBS, 1% penicillin/streptomycin, 2X L-glutamine) and transiently-transfected with a Cas9 expression vector (pLX311-Cas9) along with one of several sgRNA expression vectors (pXPR003) at a 10:1 w:w ratio using TransIT-LT1 transfection reagent (Mirus, Madison, WI, USA) (See Supplementary Table 4). 48 h after transfection, cells were trypsinized and re-plated in media containing 2.5 μM nutlin-3 (Cayman Chemical, Ann Arbor, MI, USA). Cells were expanded under nutlin-3 selection for 4 weeks to enrich for cells that had deleted WT p53. Thereafter, cells were maintained in media lacking nutlin-3. Stable p53NULL populations were derived from three independent sgRNAs, one of which was chosen for subsequent rescue experiments and genetic screens (sgTP53-4, 5′-CCCCGGACGATATTGAACAA-3′). Using the Tracking of Indels by DEcomposition algorithm[49] (TIDE - https://tide.deskgen.com), we confirmed the presence of single base insertions (57.5% of sequences, P = 0.0), and deletions of one base (10.2%, P = 3.5 x 10−59), two bases (9.8%, P = 4.4 x 10−55), or four bases (18.3%, P = 3.6 x 10−180), all of which lead to premature termination codons. Loss of endogenous p53 protein expression was determined in all populations by immunoblot. Although we could have selected a single clone for these studies, the use of cell populations decreases the likelihood of studying clone-specific effects.

Immunoblots

Cells were seeded into 6-well dishes at 2 x 105 cells/well in standard culture media. Twenty-four hours later, cells were treated with nutlin-3 (10 μM) or DMSO vehicle (0.1%) (Sigma Aldrich, St. Louis, MO, USA), and incubated for a further 24 hours. Cells were then washed in ice cold PBS (Corning, Corning, NY, USA) and protein was isolated using RIPA buffer (Cell Signaling Technology, Danvers, MA, USA) containing protease inhibitors (complete, Roche, Mannheim, Germany) and phosphatase inhibitors (PhosSTOP, Roche, Mannheim, Germany). Protein was quantified by BCA assay (Thermo Scientific, Rockford, IL, USA) and equal amounts of protein were heated to 95°C in LDS buffer (Invitrogen Carlsbad, CA, USA) containing 1% 2-mercaptoethanol (Sigma Aldrich, St. Louis, MO, USA). Protein (20 μg) was loaded onto a 4–12% Bis-Tris gel (Invitrogen, Carlsbad, CA, USA), resolved by gel electrophoresis, and transferred onto a nitrocellulose membrane using an iBlot apparatus (Invitrogen, Carlsbad, CA, USA). Membranes were blocked in PBS/0.05% Tween containing 10% milk (LabScientific, Highlands, NJ, USA) or 10% BSA (Fraction V, Fisher Scientific, Fairlawn, NJ, USA) and probed with primary antibodies targeting p53 (DO-1, Santa Cruz Biotechnology, Dallas, TX, USA), p21 (12D1, Cell Signaling Technology, Danvers, MA) or beta-actin (C4-HRP-conjugated, Santa Cruz Biotechnology, Dallas, TX, USA). HRP-conjugated secondary antibodies targeting mouse or rabbit IgG (Thermo Scientific, Rockford, IL, USA) were used in conjunction with Western Lightning Plus ECL reagent (Perkin-Elmer, Waltham, MA, USA) to visualize protein bands.

Expression of wild-type or mutant p53

A construct encoding mutant p53 p.Pro278Ala was obtained from the Dana-Farber Cancer Institute CCSB ORFeome in pDONR223. To produce a plasmid encoding WT p53, this plasmid was subjected to site-directed mutagenesis using a Quikchange Lightning II xL kit (Stratagene, Cedar Creek, TX, USA). Both p53 variants, or a Renilla luciferase control ORF were cloned into the pLX313 lentiviral destination vector via Gateway cloning (Invitrogen) and lentivirus was produced according the Broad Institute Genetic Perturbation Platform protocol (GPP; see URLs). A549 p53WT and p53NULL were then stably-infected at an MOI < 1 and selected in hygromycin B (1 mg/mL) (Invivogen, San Diego, CA, USA). p53 protein expression and activity (p21 induction) were assessed by immunoblot.

Cell viability assays

Cells were seeded at 200 cells/well in clear-bottom, opaque-walled 96-well dishes (Costar, Corning, NY, USA) in 200 μL of normal culture media and allowed to adhere overnight. Cells were then treated with nutlin-3 or etoposide (Sigma Aldrich, St. Louis, MO, USA), at 10 concentrations ranging from 20 μM to 39 nM (2-fold dilutions) and incubated for 7 days with periodic visual inspection. At the end of each experiment, media was aspirated from each well and 50 μL of a 1:1 mixture of CellTiter-Glo reagent (Promega, Madison, WI, USA) and PBS (Corning, Corning, NY, USA) was added. Plates were protected from light and incubated for 20 minutes at room temperature before being read on a Wallac EnVision plate reader (Perkin-Elmer, Waltham, MA, USA). Readings from drug-treated wells were normalized to DMSO vehicle-treated wells (set to 100% luminescence) and wells containing only media (set to 0% luminescence).

LucifeRace competition assay

A549 p53WT and p53NULL cells were stably-infected at a MOI < 1 with lentivirus encoding firefly luciferase or Renilla luciferase under control of the human EF1alpha promoter (pLX313) and selected in hygromycin B (1 mg/mL) (Invivogen). Cells were mixed at a 1:1 ratio and seeded at 2,500 cells/well in two replicate 96-well dishes in 200 μL of normal culture media and allowed to adhere overnight. Cells were then treated with compounds at 6 doses and incubated for 2 days. One plate was then subjected to a dual luciferase assay[50] and luminescence readings were obtained using a Wallac EnVision (Perkin-Elmer). Readings from drug-treated wells were normalized to DMSO vehicle-treated wells (set to 100% luminescence) and wells containing only media (set to 0% luminescence). Normalized luminescence values were then expressed as a ratio. The second replica plate was passaged at a 1/4 dilution to continue the assay. Briefly, media was aspirated, 25 μL of trypsin 0.25% (Invitrogen) was added, and plates were incubated for 5 minutes. To quench the trypsin, 175 μL of culture media was added and cells were lifted and mixed by repeated pipetting. Fifty microliter aliquots of these suspensions were then transferred to two new replicate plates, each containing 150 μL of media and fresh compound. The process of reading and re-plating was repeated every two days.

Gene synthesis and cloning

The EM7 promoter and TP53 coding sequence bearing an amber stop codon were synthesized (GenScript, Piscataway, NJ, USA) in frame with the Sh ble gene, which confers Zeocin resistance[51]. This construct was sequence confirmed and cloned into pUC57-KAN (EcoRV/HindIII) (NEB, Ipswich, MA, USA) to create the entry vector backbone template pUCTP53EV.

Oligo synthesis and tile amplification

The TP53 open reading frame was tiled using 150 base oligos each having 30 bases of complementary sequence flanking a ninety-base variable region[13]. Oligos were synthesized (CustomArray, Bothell, WA, USA) with all possible amino acid substitutions at each position. Adjacent tile oligos were manufactured on separate medium density (12K) array chips. Raw oligo tile pools were amplified by emulsion PCR (MICELLULA kit; Roboklon, Berlin, GER) using primers designed to the 30 base constant sequence regions and PfuUltra II master mix (Agilent, Santa Clara, CA, USA). Emulsions were split into seven 50 μl reactions, thermal cycled, cleaned per the MICELLULA kit instructions and then purified on a 2% E-Gel (ThermoFisher, Waltham, MA, USA).

Entry vector assembly

Primers designed to anneal to the constant region for each tile and Phusion polymerase (NEB, Ipswich, MA, USA) were used to linearize pUCTP53EV. The PCR reactions were purified on a 1% E-Gel. DpnI treated linear plasmid backbone was mixed with the relevant tile and assembled via multiplexed in vitro recombination[13] with the NEBuilder HiFi DNA Assembly Master Mix (NEB). The assembly reactions were purified, electroporated into TG1 E. coli cells (Lucigen, Middleton, WI, USA), and recovered for one hour at 37˚C in Recovery Media (Lucigen). Aliquots from the transformations were used to inoculate overnight cultures of LB containing 25 μg/mL of Zeocin (ThermoFisher, Waltham, MA, USA) to eliminate constructs harboring frameshift mutations. Cells were harvested by centrifugation and plasmid DNA was isolated using the QIAGEN Midiprep Plus kit (QIAGEN, Germantown, MD, USA). Entry libraries were verified by sequencing Nextera XT (Illumina, San Diego, CA, USA) preparations of each tile plasmid pool.

Expression library construction

The lentiviral vector pMT_BRD025 was developed by the Broad Institute Genetic Perturbation Platform (GPP). Open Reading Frames (ORF) can be cloned into this entry vector through restriction/ligation at a multiple cloning site (MCS), which allows for ORF expression under control of the human EF1α promoter. A PAC gene, which confers puromycin resistance, is driven by the SV40 promoter. To clone the TP53-MITE entry library into this vector, library DNA was digested with the restriction enzymes NheI/MluI (NEB) and ligated with pMT_BRD025 vector that was opened with these same two enzymes. To avoid bottlenecks in clone distribution, we sought to obtain 1000 bacterial colonies per variant, or 8 million colonies for the entire TP53-MITE expression library. Plasmid DNA (pDNA) was extracted from the harvested colonies using QIAGEN Maxi Prep Kits. The resulting pDNA library was sequenced via Illumina Nextera XT platform to determine the distribution of variants within the library.

Lentivirus production

Lentivirus was produced and titered by the Broad Institute GPP. The detailed protocols are available (GPP; see URLs). Briefly, 293T viral packaging cells were transfected using TransIT-LT1 transfection reagent (Mirus Bio) with the pDNA library, a packaging plasmid containing gag, pol and rev genes (psPAX2, Addgene), and an envelope plasmid containing VSV-G (pMD2.G, Addgene). Media was changed 6–8 hours after transfection and virus was harvested 24 hours thereafter. Appropriately diluted virus was used to infect A549 cells. The infected cells were then selected in puromycin and surviving cells were measured by Alamar Blue staining. Lentiviral titer was determined to be 4.4 x 105 particles/mL.

TP53-MITE library screen

A549 p53WT and p53NULL cells were infected with the TP53-MITE library using polybrene infection reagent (5 μg/mL) (Santa Cruz Biotechnology) in two independent experiments (4 x 107 cells mixed with 1.2 x 107 viral particles) and selected in puromycin (2 μg/mL) (Invivogen). Cells were then split into five equal fractions: one aliquot was saved as an early time point reference (ETP), one was treated with DMSO vehicle at 0.1% v/v, one with etoposide at 5 μM, one with nutlin-3 at 2.5 μM, and one with nutlin-3 at 5 μM. Nutlin-3- and DMSO-treated cells were trypsinized, counted, and re-seeded every 3 days for 12 days, whereas etoposide-treated cells were harvested and counted only on Day 12. Aliquots comprising approximately 8 x 106 cells were frozen in PBS (Corning) at each harvest. The following samples were subjected to gDNA isolation and subsequent analysis: p53WT ETP, p53WT + nutlin-3 (2.5 μM), p53NULL ETP, p53NULL + nutlin-3 (5 μM), p53NULL + etoposide (5 μM). gDNA was isolated using QIAamp DNA Blood Midi kits (QIAGEN). Spectrophotometric analysis of the gDNA indicated that between 2.2 x 106 and 5.5 x 106 genome equivalents were recovered from each condition.

ORF purification from gDNA

Twelve PCR reactions were performed for each gDNA sample. The volume of each PCR reaction was 100 μL and contained ~3 μg of gDNA. Herculase II (Agilent Genomics) was used as the DNA polymerase. All 12 PCR reactions for each gDNA sample were pooled, concentrated with a PCR cleanup kit (QIAGEN), loaded onto a 1% agarose gel, and separated by gel electrophoresis. Bands of the expected size were excised, and DNA was purified first using a QIAquick kit (QIAGEN) then an AMPure XP kit (Beckman Coulter).

Nextera sequencing

Reactions were performed according to the Illumina Nextera XT protocol. For each sample, we set up 6 Nextera reactions, each with 1 ng of purified ORF DNA. Each reaction was indexed with unique i7/i5 index pairs. After the limited-cycle PCR step, the Nextera reactions were purified with AMPure XP kit. All samples were then pooled and sequenced using an Illumina Hiseq4000 platform using two reads, each 150 bases in length.

TP53-MITE library screen data analysis

Hiseq4000 data were processed with ORFcall software developed by the Broad Institute[15] (See Supplementary Note) and aligned to the TP53 reference sequence. The number of reads corresponding to each variant were then tallied. At each codon position, there were counts for all 20 amino acids and a stop codon. At 378 codon positions there were counts for synonymous codons (silent mutations) (Supplementary Table 2). The raw read counts were normalized to the fraction of counts at each codon position. Experimental replicates yielded highly correlated read counts for each allele (Supplementary Fig. 4a, d, g). To assess the level of enrichment or depletion of each amino acid variant, we calculated the log2-fold change in fractional read counts relative to the early time point samples for all non-wild-type residues and synonymous variants, then we averaged and standardized these values to produce Z-scores (Supplementary Table 3, Fig. 2, and Supplementary Fig. 4 b–c, e–f, h–i). In the attached data files, synonymous variants are represented by the amino acid letter “B” and stop codons are represented by the letter “Z”. A Combined Phenotype Score was calculated as: (Z-score(p53WT + nutlin) + Z-score(p53NULL + nutlin) − Z-score(p53NULL + etoposide))/3 ± SD(+ Z-score(p53WT + nutlin), + Z-score(p53NULL + nutlin), − Z-score(p53NULL + etoposide))/√ 3.

Assignment of mutation probability based on trinucleotide mutation context

Each base in the wild-type TP53 cDNA sequence was systematically changed to every other base in silico, and the resulting sequences were translated using the standard genetic code. By considering the nucleotide before and after the altered base, each translated mutant p53 protein was assigned to its trinucleotide mutation type. The 192 trinucleotide mutation types were collapsed to 96 types by considering base changes from the perspective of the complementary strand read in the opposite direction. For example, 5′-A[G>T]C-3′ mutations were considered equivalent to 5′-G[C>A]T-3′ mutations. Mutation probabilities were then assigned to each variant using all 30 signatures in the COSMIC mutational signatures database, as well as a transcription strand specific version of Signature 4 (Signature 4*, Supplementary Table 5). Mutation probabilities for nucleotide changes that yield the same amino acid change were then added together as these events were assumed to be rare and independent in a population.

Systematic identification of tissue-of-origin-enriched TP53 mutations

Enrichment analyses were performed using the IARC R18 database (IARC; see URLs)[2,32]. Fisher’s exact tests were performed on each mutant allele in each tissue type that was represented by greater than 300 samples, and a Bonferroni correction was applied to correct for multiple-hypothesis testing. Enrichments were considered significant at an α level of 0.05. These mutation type-tumor type relationships were then validated in an independent dataset (AACR GENIE; see URLs)[35] using a nominal P value cutoff of 0.05.

Modeling the TP53 mutational spectrum

We created a generalized linear model in R (3.3.3 GUI 1.69 Mavericks build) using the glm function, assuming a Poisson noise distribution and a logarithmic link function. This function was used to predict the mutation count for each TP53 allele in the IARC database using inputs for all three phenotypic screens and projections of mutational Signatures 1 and 5, which are associated with aging and have been found in all tumor types, as well as Signatures 2 (APOBEC), 4* (smoking), 6 (MMR deficiency), 7 (UV), 13 (APOBEC), and 24 (aflatoxin), which were identified as contributing tissue-of-origin selective TP53 mutations.

To derive the models, the following calls were used:

Where SigX represents the mutation probability for each allele under the influence of Signature X, as defined in the COSMIC mutational signatures database, and Z_DN, Z_LOF, and Z_WT represent the Z-score for each allele in the p53WT + nutlin-3, p53NULL + nutlin-3, and p53NULL + etoposide screens, respectively. As outlined in Supplementary Fig. 6, each of these parameters significantly contributed to the combined model (P values range from 4.74 x 10−5 to < 2.2 x 10−16), and all parameters except Z_WT produced positive coefficients, which was expected, given that LOF alleles were depleted in the etoposide screen. Model inputs and predictions for each amino acid variant are reported in Supplementary Table 3.

In addition to testing the predictive accuracy of the combined model on a separate database of mutations than the one used for training (GENIE)[35], we also controlled for the possibility of the models being over-fit to correlated noise in the two mutation databases by using 10-fold cross-validation. Briefly, we created 10 complementary data partitions (folds) using the mutation counts from the IARC R18 database. For each fold, we used 90% of the alleles to train a mutational signatures model, a phenotypic selection model, and a combined model. We then used the outputs of these models to predict the mutation counts for the remaining 10% of alleles we had set aside as a test set. We found that the contributions made by each of the input variables to the combined model were similar across all folds, as well as in the final model, and that the descriptive accuracy was similar to the predictive accuracy, suggesting that the coefficients had not been over-fit (Supplementary Fig. 6).