Author response:
The following is the authors’ response to the original reviews.
Public Reviews:
Reviewer #1 (Public Review):
Summary:
Kimura et al performed a saturation mutagenesis study of CDKN2A to assess the functionality of all possible missense variants and compare them to previously identified pathogenic variants. They also compared their assay result with those from in silico predictors.
Strengths:
CDKN2A is an important gene that modulates cell cycle and apoptosis, therefore it is critical to accurately assess the functionality of missense variants. Overall, the paper reads well and touches upon major discoveries in a logical manner.
Weaknesses:
The paper lacks proper details for experiments and basic data, leaving the results less convincing. Analyses are superficial and do not provide variant-level resolution.
We thank the reviewer for their comments. We have updated the manuscript to include additional detail of experimental methods and variant level resolution of data and analyses. We have also conducted additional analyses to compare variant classifications using a gamma generalized linear model and log2 normalized fold change, establish the effect of low variant coverage on variant functional classifications, determine the performance of combining multiple in silico predictions, and determine the prevalence of functionally deleterious variants in gnomAD and functionally deleterious variants of uncertain significance in ClinVar compared all CDKN2A missense variants.
Reviewer #2 (Public Review):
This study describes a deep mutational scan across CDKN2A using suppression of cell proliferation in pancreatic adenocarcinoma cells as a readout for CDKN2A function. The results are also compared to in silico variant predictors currently utilized by the current diagnostic frameworks to gauge these predictors' performance. The authors also functionally classify CDKN2A somatic mutations in cancers across different tissues.
This study is a potentially important contribution to the field of cancer variant interpretation for CDKN2A, but is almost impossible to review because of the severe lack of details regarding the methods and incompleteness of the data provided with the paper. We do believe that the cell proliferation suppression assay is robust and works, but when it comes to the screening of the library of CDKN2A variants the lack of primary data and experimental detail prevents assessment of the scientific merit and experimental rigor.
We are grateful for the opportunity to clarify our experimental methods and to provide additional data in the revised manuscript. The manuscript has been updated to include, among other changes, additional information on assay design, analysis of variant representation in the library, inclusion of primary data with variant level resolution, and a comparison of variant classifications using a gamma generalized linear model and log2 normalized fold change.
Recommendations for the authors:
Reviewer #1 (Recommendations For The Authors):
Major issues:
(1) Can the pathogenicity values of individual amino acid changes be opened to the public? It would serve as a valuable asset to the community.
Thank you for your suggestion. We are happy to provide this information. Individual variant data and functional classifications from the functional assay are given in Appendix 1-table 4.
(2) In the method section, it is not clear (at least to the reviewer) whether the protocol describing the construction of the CDKN2A missense library was provided.
Thank you for your comment. We have included additional information in the manuscript describing construction of the CDKN2A missense library.
“CDKN2A expression plasmid libraries
Codon-optimized CDKN2A cDNA using p16INK4A amino acid sequence (NP_000068.1), was designed (Appendix 1-table 12) and pLJM1 containing codon optimized CDKN2A (pLJM1-CDKN2A) generated by Twist Bioscience (South San Francisco, CA). 156 plasmid libraries were then synthesized by using pLJM1-CDKN2A, such that each library contained all possible 20 amino acids variants (19 missense and 1 synonymous) at a given position, generating 500 ng of each plasmid library (Twist Bioscience, South San Francisco, CA). The proportion of variant in each library was shown in Appendix 1-table 2. Variants with a representation of less than 1% in a plasmid library were individually generated using the Q5 Site-Directed Mutagenesis kit (New England Biolabs, Ipswich, MA; catalog no. E0552), and added to each library to a calculated proportion of 5%. Primers used for site-directed mutagenesis are given in Appendix 1-table 13. Each library was then amplified to generate at least 5 ug of plasmid DNA using QIAGEN Plasmid Midi Kit (QIAGEN, Germantown, MD; catalog no. 12143).”
(3) The paper lacks basic experimental results. The results cover almost all possible missense variants, but it would be clearer if actual coverage values used for calculating relative enrichment were shown. Are all variants well covered? Isn't there any spurious signal due to low coverage? How many times were the experiments performed? Also, how many cells were used, what was the expected MOI, and what proportion of harvested cells is thought to have a single variant? How can you distinguish the effect of a single variant from a multiple variants effect?
We thank the reviewer for their comment. We have provided additional information in the manuscript to address these issues. Briefly, in response to each issue:
(1) We have provided read count data for all variants, used to determine functional classifications based on either gamma generalized linear model or normalized fold change, in Appendix 1-table 4.
(2) To assess if low variant coverage resulted in spurious signals, we compared prevalence of functionally deleterious classifications among variants binned by coverage in the Day 9 cell pool. We did not identify any statistically significant differences based on variant coverage.
“We also determined whether underrepresentation in the cell pool at Day 9 affected variant functional classifications. Fifty-three of 2,964 missense variants (1.8%) were present in the cell pool at Day 9 of the first assay replicate (experiment 1) at < 2%, as determined by the number of sequence reads supporting the variant (Figure 2 -figure supplement 4A, Appendix 1-table 4). There was no statistically significant difference in the proportion of variants classified as functionally deleterious for variants present in less than 2% of the cell pool at Day 9 (12 of 53 variants; 22.6%), and variants present in more than 2% of the cell pool (496 of 2,911 variants; 17.0%) (P value = 0.28) (Figure 2 -figure supplement 4B). We also found no significant differences in the proportion of variants classified as functionally deleterious for variants present in more than 2% of the cell pool at Day 9 when variants were binned in 1% intervals (Figure 2 -figure supplement 4B).”
(3) The assay was repeated in duplicate for 28 CDKN2A residues. For the remaining 128 residues of CDKN2A, the assay was completed once. We found good agreement between variant classifications in assay repeats. We have added to the text as follows:
“To confirm the reproducibility of our variant classifications, 28 amino acid residues were assayed in duplicate, and variants classified using the gamma GLM. The majority of missense variants, 452 of 560 (80.7%), had the same functional classification in each of the two replicates (Figure 2 -figure supplement 3A and B, Appendix 1-table 4).”
We have also added discussion of this study limitation to the manuscript:
“We repeated our functional assay twice for 28 CDKN2A residues. For the remaining 128 residues of CDKN2A, the functional assay was completed once. While we found general agreement between functional classifications from each replicate for the 28 residues assayed in duplicate, additional repeats for each residue are necessary to determine variability in variant functional classifications.”
(4) We have added additional information about the number of cells used for transduction and MOI to the method section:
“Lentiviral transduction
PANC-1 cells were used for CDKN2A plasmid library and single variant CDKN2A expression plasmid transductions. PANC-1 cells previously transduced with pLJM1-CDKN2A (PANC-1CDKN2A) and selected with puromycin were used for CellTag library transductions. Briefly, 1 x 105 cells were cultured in media supplemented with 10 ug/ml polybrene and transduced with 4 x 107 transducing units per mL of lentivirus particles. Cells were then centrifuged at 1,200 x g for 1 hour. After 48 hours of culture at 37oC and 5% CO2, transduced cells were selected using 3 µg/ml puromycin (CDKN2A plasmid libraries and single variant CDKN2A expression plasmids) or 5 µg/ml blasticidin (CellTag plasmid library) for 7 days. Expected MOI was one. After selection, cells were trypsinized and 5 x 105 cells were seeded into T150 flasks. DNA was collected from remaining cells and this sample was named as (Day 9). T150 flasks were cultured until confluent and then DNA was collected. The time for cells to become confluent varied for each amino acid residue (Day 16 – 40, Appendix 1-table 5).”
(5) Our assay was not designed to distinguish multiple variant effects. However, we do not anticipate multiple transductions to significantly impact variant classifications in our assay. We found that our functional classifications were consistent with previously reported classifications:
“In general, our results were consistent with previously reported classifications. Of variants identified in patients with cancer and previously reported to be functionally deleterious in published literature and/or reported in ClinVar as pathogenic or likely pathogenic (benchmark pathogenic variants), 27 of 32 (84.4%) were functionally deleterious in our assay (Figure 2B, Figure 2 -figure supplement 1B and 1C, Appendix 1-table 4) (Chaffee et al., 2018; Chang et al., 2016; Horn et al., 2021; Hu et al., 2018; Kimura et al., 2022; McWilliams et al., 2018; Roberts et al., 2016; Zhen et al., 2015). Five benchmark pathogenic variants were characterized as indeterminate function, with log2 P values from -19.3 to -33.2. Of 156 synonymous variants and six missense variants previously reported to be functionally neutral in published literature and/or reported in ClinVar as benign or likely benign (benchmark benign variants), all were characterized as functionally neutral in our assay (Figure 2B, Figure 2 -figure supplement 1B and 1C, Appendix 1-table 4) (Kimura et al., 2022; McWilliams et al., 2018; Roberts et al., 2016). Of 31 VUSs previously reported to be functionally deleterious, 28 (90.3%) were functionally deleterious and 3 (9.7%) were of indeterminate function in our assay. Similarly, of 18 VUSs previously reported to be functionally neutral, 16 (88.9%) were functionally neutral and 2 (11.1%) were of indeterminate function in our assay, (Figure 2B, Figure 2 -figure supplement 1B and 1C, Appendix 1-table 4).”
(4) Comparison of functional classifications (shown in Figure 3) from this study and other in silico tools is superficial. The analysis is based on the presumption that their result is gold-standard, thereby calculating the sensitivity, accuracy, and PPV of individual predictors. But apparently, this won't be true, so it would be more reasonable to check the "correlation" of the study results and other predictors: e.g. which variants show consistent results between this study and other predictors? Are there any indicators of consistent vs inconsistent results? How does the consistency change by protein sequences or domains? Etc
Thank you for your comment. We have added additional analysis to our manuscript comparing our functional classifications with in silico variant effect predictions. Specifically, we have included analysis combining multiple predictors:
“We also tested the effect of combining multiple in silico predictors. 904 missense variants had in silico predictions from all 7 algorithms. The remaining 2,060 missense variants had in silico predictions from 5 algorithms. Of variants with in silico predictions from all 7 algorithms, 378 (41.8%) had predictions of deleterious or pathogenic effect from a majority of algorithms (≥ 4), and of these, 137 (36.2%) were functionally deleterious in our assay. Similarly, of 2,060 missense variants that had in silico predictions from 5 algorithms, 1107 (53.7%) had predictions of deleterious or pathogenic effect from a majority of algorithms (≥ 3), of which, 361 (32.6%) were functionally deleterious in our assay (Appendix 1-table 7).”
(5) Similarly, Figure 4 does not deliver much information, either. Rather than delivering a simple summary, it would be more informative if deeper analyses were conducted. e.g., do pathogenic variants show higher frequency among patients, or higher variant frequency in tumors (if data were available).
We have included additional analysis of somatic alterations in the manuscript. We found pathogenic/likely pathogenic somatic mutations were enriched in patients. This was also the case for somatic mutations that were classified as functionally deleterious in our assay. We also found statistically significant depletion of functionally deleterious mutations in colorectal adenocarcinoma. Interestingly, no patients with a somatic mutation in a mismatch repair gene had a functionally deleterious CDKN2A missense somatic mutation. However, this observation was not statistically significant. Future studies will determine whether CDKN2A and MMR gene somatic mutations are mutually exclusive in colorectal adenocarcinoma.
“We found that 34.2% - 53.4% of unique missense somatic mutations classified as functionally deleterious, with 61.4% - 67.6% of patients having a functionally deleterious somatic mutation (Figure 4A, Appendix 1-table 9). As with functionally deleterious variants, functionally deleterious missense somatic mutations were also not distributed evenly across CDKN2A, being enriched within the ankyrin repeat 3 (Figure 4B, Appendix 1-table 9). We found that 32.4% - 50.0% of all functionally deleterious missense somatic mutations occurred within ankyrin repeat 3, with 48.0% - 58.0% of patients in each cohort having a functionally deleterious missense somatic mutation in this domain. Notably, 65.7% - 76.0% of functionally deleterious missense somatic mutations in this domain were in residues 80-89 (Appendix 1-table 9).”
“We were also able to determine the functional classification of CDKN2A missense somatic mutations in COSMIC, TCGA, JHU, and MSK-IMAPCT by cancer type. We found that 22.2% - 100% of CDKN2A missense somatic mutations were functionally deleterious depending on cancer type (Figure 4-figure supplement 2A-D). When considering missense somatic mutation reported in any database, there was a statistically significant depletion of functionally deleterious mutations in colorectal adenocarcinoma (20.4%; adjusted P value = 5.4 x 10-9) (Figure 4C). As the proportion of missense somatic mutations that were functionally deleterious was less in colorectal carcinoma compared to other types of cancer, we assessed whether somatic mutations in mismatch repair genes (MLH1, MLH3, MSH2, MSH6, PMS1, and PMS2) were associated with the functional status of CDKN2A missense somatic mutations. Thirty-five patients in COSMIC had a CDKN2A missense somatic mutation, of which 12 (34.3%) had a somatic mutation in a mismatch repair gene. We found that no patients with a somatic mutation in a mismatch repair gene had a functionally deleterious CDKN2A missense somatic mutation compared to 6 of 23 samples (26.1%) without a somatic mutation in a mismatch repair gene (P value = 0.062).”
(6) It would be helpful to validate the neutral variants set. Are variants of UK biobank or gnomAD enriched on neutral population? Are synonymous variants exclusively found in neutral populations?
Thank you for the suggestion. All synonymous variants were found to functionally neutral in our assay. We also assessed VUSs from gnomAD and found a lower prevalence of functionally deleterious variants compared to all CDKN2A variants and CDKN2A missense somatic mutations:
“The Genome Aggregation Database (gnomAD) v4.1.0 reports 287 missense variants in CDKN2A, including the 13 pathogenic, 4 likely pathogenic, 3 likely benign, 3 benign, and 264 VUSs classified using ACMG variant interpretation guidelines (Figure 5A, Figure 5B, and Appendix 1-table 10). Of the 264 missense VUSs, 177 were functionally neutral (67.0%), 56 (21.2%) were indeterminate function, and 31 (11.7%) were functionally deleterious in our assay using the gamma GLM for classification (Figure 5C).”
(7) They used a pancreatic cancer cell line and assayed for cell proliferation. The limitations of this method and the possibility of complementing the limitations should be discussed.
Thank you for the suggestion. We have added discussion of this limitation to our manuscript:
“We characterized variants based upon a broad cellular phenotype, cell proliferation, in a single PDAC cell line. It is possible that CDKN2A variant functional classifications are cell-specific and assay-specific. Our assay may not encompass all cellular functions of CDKN2A and an alternative assay of a specific CDKN2A function, such as CDK4 binding, may result in different variant functional classifications. Furthermore, CDKN2A variants may have different effects if alternative cell lines are used for the functional assay. However, cell-specific effects appear to be limited. In our previous study, we characterized 29 CDKN2A VUSs in three PDAC cell lines, using cell proliferation and cell cycle assays, and found agreement between all functional classifications (Kimura et al., 2022).”
Minor issues:
(1) Figures 2B, C: it would be more intuitive to plot significance by logging p-values than raw p-values.
We used log2 P value (or log2 normalized fold change) for figures in the manuscript as appropriate.
(2) Figure 2D: annotate protein domain information at the side. Supplementary Figure 2 shows the domains but it would be more informative to show it in Figure 2D heatmap.
Thank you for the suggestion, we have annotated protein domain information on the left side of the heatmap in (the now) Figure 2C.
Reviewer #2 (Recommendations For The Authors):
Major Concerns:
(1) How many replicates of the screen were performed? It seems like only one library infection/ proliferation assay was done. If so this is insufficient to obtain any idea of the uncertainty of measurement for each variant.
The assay was repeated in duplicate for 28 CDKN2A residues. For the remaining 128 residues of CDKN2A, the assay was completed once. We found good agreement between variant classifications in assay repeats. We have added to the text as follows:
“To confirm the reproducibility of our variant classifications, 28 amino acid residues were assayed in duplicate, and variants classified using the gamma GLM. The majority of missense variants, 452 of 560 (80.7%), had the same functional classification in each of the two replicates (Figure 2 -figure supplement 3A and B, Appendix 1-table 4).”
We have also added discussion of this study limitation to the manuscript:
“We repeated our functional assay twice for 28 CDKN2A residues. For the remaining 128 residues of CDKN2A, the functional assay was completed once. While we found general agreement between functional classifications from each replicate for the 28 residues assayed in duplicate, additional repeats for each residue are necessary to determine variability in variant functional classifications.”
(2) The count data from the experiment and NGS pipeline to call variants need to be provided for each replication (i.e. the counts that were fed into the gamma model)
Accompanying this should be information about the depth of sequencing of the cells, the number of cells infected with the library, and standard metrics for pooled screens.
Quality metrics regarding the representation and completeness of the TWIST library need to be provided. See Brenan et al. Cell Reports (2016) Supplemental Figure 1
Thank you for your suggestion. We are happy to provide this additional information. Sequence read counts for each variant are given in Appendix 1-table 4. We have provided addition detail in the methods section on functional assay, including number of cells infected with each library:
“Lentiviral transduction
PANC-1 cells were used for CDKN2A plasmid library and single variant CDKN2A expression plasmid transductions. PANC-1 cells previously transduced with pLJM1-CDKN2A (PANC-1CDKN2A) and selected with puromycin were used for CellTag library transductions. Briefly, 1 x 105 cells were cultured in media supplemented with 10 ug/ml polybrene and transduced with 4 x 107 transducing units per mL of lentivirus particles. Cells were then centrifuged at 1,200 x g for 1 hour. After 48 hours of culture at 37oC and 5% CO2, transduced cells were selected using 3 µg/ml puromycin (CDKN2A plasmid libraries and single variant CDKN2A expression plasmids) or 5 µg/ml blasticidin (CellTag plasmid library) for 7 days. Expected MOI was one. After selection, cells were trypsinized and 5 x 105 cells were seeded into T150 flasks. DNA was collected from remaining cells and this sample was named as (Day 9). T150 flasks were cultured until confluent and then DNA was collected. The time for cells to become confluent varied for each amino acid residue (Day 16 – 40, Appendix 1-table 5). DNA was extracted from PANC-1 cells using the PureLink Genomic DNA Mini Kit (Invitrogen, Carlsbad, CA; catalog no. K1820-01). The assay for CellTag library was repeated in triplicate. We repeated our CDKN2A assay in duplicate for 28 residues. For the remaining 128 CDKN2A residues the assay was completed once.”
We have also provided additional information on the TWIST library:
“CDKN2A expression plasmid libraries
Codon-optimized CDKN2A cDNA using p16INK4A amino acid sequence (NP_000068.1), was designed (Appendix 1-table 12) and pLJM1 containing codon optimized CDKN2A (pLJM1-CDKN2A) generated by Twist Bioscience (South San Francisco, CA). 156 plasmid libraries were then synthesized by using pLJM1-CDKN2A, such that each library contained all possible 20 amino acids variants (19 missense and 1 synonymous) at a given position, generating 500 ng of each plasmid library (Twist Bioscience, South San Francisco, CA). The proportion of variant in each library was shown in Appendix 1-table 2. Variants with a representation of less than 1% in a plasmid library were individually generated using the Q5 Site-Directed Mutagenesis kit (New England Biolabs, Ipswich, MA; catalog no. E0552), and added to each library to a calculated proportion of 5%. Primers used for site-directed mutagenesis are given in Appendix 1-table 13. Each library was then amplified to generate at least 5 ug of plasmid DNA using QIAGEN Plasmid Midi Kit (QIAGEN, Germantown, MD; catalog no. 12143).”
(3) It is unclear when barcode abundance is assessed in the cell proliferation assay/in the screen. The exact timepoints of "before and after in vitro culture" (line 91) need to be clarified in the text.
We are happy to clarify. We collected DNA on Day 9 post transfection and at confluency. Day of confluency for each residue is detailed in Appendix 1-table 5. The text of the manuscript has been updated appropriately.
(4) Is "before" day 9, as detailed in Figure 1 source data 1? If so, it is misleading to state that the experiment is in culture for 14 days but call day 9 "before... in vitro culture."
The "before" sample should be obtained immediately after viral infection and selection with the library to provide a representation of library representation.
We apologize for your confusion. We have clarified in the text and figures that our baseline measurement was at Day 9 post transfection. We also determined whether the proportion of each variant is maintained in the Day 9 cell pool compared to the amplified plasmid library for three CDKN2A amino acid residues (p.R24, p.H66, and p.A127) and updated the manuscript text:
“To confirm that the representation of each variant was maintained after transduction, we transduced three lentiviral libraries (amino acid residues p.R24, p.H66, and p.A127) individually into PANC-1 cells and determined the proportion of each variant in the amplified plasmid library and in the cell pool at Day 9 post-transduction. The proportion of each variant in the amplified plasmid library and in the cell pool at Day 9 were highly correlated (Figure 1 -figure supplement 2C and D, Appendix 1-table 3).”
(5) There is no information regarding the function of each variant, aside from just a p-value resulting from the final analysis with the gamma model. Some variants may cause loss of function, others may be neutral while others may be gain of function. Simply providing a p-value is not sufficient. The standard in the field is to provide a function score/ test-statistic giving the sign and magnitude of the effect. For proliferation assays at least a ratio of fold-change of (mut/ synonymous)[day 14] vs (mut/synonymous)[baseline] should be provided.
Thank you for your comment. We have provided read counts, P values, and functional classifications for each variant using the gamma GLM in Appendix 1-table 4. We have also analyzed variants using log2 normalized fold change. This data is presented in the text and compared to our classifications with the gamma GLM. We have provided normalized fold change and resulting classification for each variant in Appendix 1-table 6.
(6) A plot of the distribution of function scores for all variants is needed. This will serve as an effective visual to distinguish the control variants from those that are functionally deleterious or benign/neutral (see Findlay et al. Nature (2018) Figure 3A for an example visual).
Thank you for your suggestion. We have provided additional figures to visualize distribution of assay outputs using the gamma GLM in Figure 2 -figure supplement 1.
(7) Synonymous variants are used as a proxy for WT per variant library, but do all the synonymous variants truly behave like WT CDKN2A in their ability to suppress cell proliferation? A plot of the distribution of synonymous variant function relative to WT CDKN2A function would be effective here.
All 156 synonymous variants suppressed cell proliferation and were classified as functionally neutral in our assay using the gamma GLM. The manuscript has been updated to reflect this:
“Of 156 synonymous variants and six missense variants previously reported to be functionally neutral in published literature and/or reported in ClinVar as benign or likely benign (benchmark benign variants), all were characterized as functionally neutral in our assay (Figure 2B, Figure 2 -figure supplement 1B and 1C, Appendix 1-table 4)”
(8) The gamma generalized linear model is not commonly used to analyze the results of saturation mutagenesis screens. Please provide a justification for the use of this analysis method vs using log fold change as other dms scan studies have done (PMID: 27760319, PMID: 30224644).
Thank you for this important suggestion. We are happy to provide additional information. We used a gamma GLM to functionally characterize CDKN2A variants as it does not rely on an annotated set of pathogenic and benign variants to determine classification thresholds. Instead, classification thresholds are determined using the change in representation of 20 non-functional barcodes in a pool of PANC-1 cells stably expressing CDKN2A after a period of in vitro growth. As a gamma GLM is not commonly used for saturation mutagenesis screens, as noted by the reviewer, we also classified variants using log2 normalized fold change. We compared variant functional classifications using the gamma GLM and log2 normalized fold change and in general we found agreement between both methods with 98.5% of missense variants classified as functionally deleterious using a gamma GLM, similarly classified using log2 normalized fold change. We have updated the text to reflect this reasoning and additional analysis.
(9) The statistical methods used to calculate enrichment of deleterious variants per region of CDKN2A (Figure 2 supplement 1B; lines 163-168) are not described anywhere in the paper. Additionally, the same statistical analysis is not applied to the variants in the subregions near the ankyrin repeats (lines 168-172).
We are happy to clarify and have added text to the methods section:
“Z-tests with multiple test correction performed with the Bonferroni method was used in the following comparisons: 1) proportion of functionally deleterious variants present in < 2% of the cell pool and ≥ 2% of the cell pool at Day 9 binned in 1% intervals, 2) proportion of variants in each domain predicted to have deleterious or pathogenic effect by the majority of algorithms, 3) proportion of functionally deleterious variants in each domain, and 4) proportion of functionally deleterious missense variants and somatic mutations.”
Minor:
(1) Please review the manuscript for spelling and grammatical errors.
Sure.
