Prime editor-based high-throughput screening reveals functional synonymous mutations in human cells

Authors

1. Xuran Niu
2. Wei Tang
3. Yongshuo Liu
4. Binrui Mo
5. Ying Yu
6. Ying Liu
7. Wensheng Wei

Extended Data Fig. 1 Design of the epegRNA library for screening functional synonymous mutations and development of the experimental system using PEmax.

a, Diagram of the epegRNA^eBAR library structure. RTT: reverse transcription template, PBS: primer binding site. b, Principles and workflow for designing epegRNAs within the library. Searches for NGG PAMs are conducted on the representative transcript of each gene to determine the spacer sequence, PBS sequence, and RTT, followed by a first filtering step. Comprehensive searches along the coding region are then performed to achieve saturation mutagenesis of synonymous mutations, with linker sequences subsequently designed and a secondary filter applied, resulting in the final epegRNA sequences. c, Distribution of expression levels and essentiality of genes within the library in HCT116 cells. Red dots represent the 11 genes targeted for complete saturation mutagenesis, blue dots represent the 56 genes targeted for synonymous saturation mutagenesis, and gray dots represent all other genes. d, Histogram of the number of epegRNA designed for mutations in the library. e, Statistics of the types of amino acid substitutions for the 11 genes targeted for saturation mutagenesis. The y-axis represents the original amino acid, the x-axis represents the amino acid post-PE editing, and the values represent the number of epegRNAs designed for each substitution. The histogram on the right illustrates the coverage of the corresponding amino acids, indicating the proportion of amino acids that can be targeted among all amino acids in the 11 genes. f, Clustered heatmap displaying Pearson correlation coefficients of the whole transcriptome between HCT116-PEmax cells and wild-type (WT) HCT116 cells, as well as after the addition of nontargeting (NT) epegRNA and AAVS1-targeting epegRNA. g, Graph showing the changes in editing efficiency over a 28-day period of continuous culture in HCT116-PEmax cells with the addition of two test epegRNAs targeting different sites. The data is presented as the mean ± s.d. (n = 3 biological replicates).

Extended Data Fig. 2 Decoding and ROC analysis of the epegRNA^eBAR library.

a, Diagram illustrating the genomic PCR process performed on the cell library post-screening. b, Algorithm description of ZFC-eBAR (see Methods for details). c, The x-axis represents the number of genes corresponding to the positive controls (epegRNAs that introduces nonsense and frameshift mutations on essential genes), while the y-axis represents the AUC calculated using these positive controls as well as negative controls (AAVS1-targeting and nontargeting epegRNA). The score used to calculate AUC is the screen score. d, Similar to c, but nonsense/frameshift mutations at the same site are regarded as the same mutation. The score used to calculate AUC is the Top score. e, ROC analysis of epegRNAs in highly essential genes (Chronos score < −2.5) based on different scores. The red point represents the threshold we selected. The red line represents the ROC curve for the top1 score of epegRNA, the dark blue line represents the ROC curve for the top2 mean score of epegRNA, and the green line represents the ROC curve for the all mean score of epegRNA. f, Spearman correlation of zLFC for the top 2 ranked epegRNAs in highly essential genes (Chronos score < −2.5).

Extended Data Fig. 3 Evaluation of screening reproducibility and detection sensitivity.

a, Correlation of zLFC values between two eBARs replicates for epegRNA with the third eBAR’s |LFC| ≥ 0. b, Correlation of zLFC values between two eBARs for epegRNAs with |AGCA LFC| ≥ 1 and |CACT LFC| ≥ 1. c, Pearson correlation analysis of LFC values between two eBARs across varying thresholds of absolute LFC values from the third eBAR. The x-axis indicates the absolute LFC threshold applied to each eBAR, while the y-axis represents the Pearson correlation coefficient between the remaining two eBARs. d, Correlation of LFC values between two eBARs for epegRNAs with the third eBAR’s |LFC| ≥ 1. e, Correlation of LFC values between two eBARs for epegRNAs with the third eBAR’s |LFC| ≥ 0. For panels a, b, d, and e, Pearson correlation coefficient (r) is labeled. N indicates the number of epegRNAs shown in these figures (for panels a and e, excluding those with zero sequencing read counts and black lines indicate density-based contours). f–g, Correlation analysis of enriched mutations in the HCT116 screen across three eBARs, based on LFC values (f) or zLFC values (g). Diagonal panels show histograms of value distributions, lower left panels show scatter plots of pairwise correlations, and upper right panels present Pearson correlation coefficient (r). The dataset includes 1,717 epegRNAs with nonsynonymous mutations and 417 epegRNAs with synonymous mutations. h, Boxplot of the relationship between the degree of reduction in counts of simulated deleterious variants in the experimental group (x-axis) and the calculated |screen score| (y-axis). n = 492 per group (number of epegRNAs). The gray dashed lines correspond to thresholds of 3 and 2.2. Boxplots are depicted as follows: the center line represents the median, the box limits denote the upper and lower quartiles, and whiskers extend to 1.5 times the interquartile range. i, Corresponding to h, the x-axis represents the average counts of these simulated deleterious variants in day 0, and the y-axis indicates the percentage reduction in day 35. Green labels highlight points where the |screen score| after reduction exceeds the selection threshold (3).

Extended Data Fig. 4 Performance of synonymous mutations in human homologous genes of the yeast.

a, Comparison of fitness score distributions for synonymous mutations across three groups: 67 saturation-mutated genes (synonymous mutations only), 11 saturation-mutated genes (including all mutation types), and human homologous genes of the yeast (excluding RPL39 due to limited data). P values were calculated using a two-sided Wilcoxon test. b, Distribution of essentiality for human homologous genes of yeast in HCT116 cells, data of RPL29 are missing in DepMap. The corresponding positions are labeled in purple. c, Shown in order as in b, distribution of fitness scores for synonymous, nonsense and frameshift mutations for each gene. From left to right, the mutation counts for each group are 46/5 (RPS7), 269/6 (TSR2), 459/6 (BUD23), 167/6 (ATP6V1F), 791/6 (PAF1), 233/6 (GET1), 131/6 (VMA21), 162/6 (LSM1), 193/6 (UBE2B), 323/6 (INO80C), 487/6 (PRPS2) and 80/6 (RPL29), with n = 494 for AAVS1 and n = 981 for Nontargeting. Boxplots are depicted as follows: the center line represents the median, the box limits denote the upper and lower quartiles, and whiskers extend to 1.5 times the interquartile range.

Extended Data Fig. 5 Comprehensive analysis of different mutation types and their effects in saturation mutagenesis.

a, Correlation of BLOSUM62 scores with the screening results for 11 genes. Mean scores were compared to nontargeting controls, with mutations showing adjusted P < 0.05 (two-tailed Student’s t test with Benjamini–Hochberg correction) marked in black. Only mutations with counts ≥10 were included. b, Upper panel displays a heatmap of screen scores for all mutations in RAN. Light gray indicates amino acids not included in the library. Synonymous mutations are marked with dots, and dashed lines denote domain boundaries. The lower panel illustrates domain annotations for RAN. c, Top scores for missense and synonymous mutations are compared across all regions of RAN, including G3 Switch II and Q69, n indicates the number of mutations. d, Comparison of top screen scores from our study with scores from ref. ²². Spearman correlation coefficient = 0.3076786. e, Detailed Mutation Analysis in POLR2L. Middle panel: Heatmap of screening scores for all mutations in POLR2L. Light gray indicates amino acids not designed in the screening library. Dots indicate synonymous mutations. Dashed lines represent four Zn²⁺ binding sites. Lower panel: Average top scores across each position, with three designed Zn²⁺ binding sites highlighted in red. Left panel: Average top scores by amino acid, with the two most depleted ones shown in orange. f, Cell fitness effects of different amino acid (aa) substitutions (statistical methods as a). g, The distributions of fitness scores for mutation types with a significant impact on cell fitness. From left to right, the mutation counts are 1,475, 182, 258, 113, 233, 230, and 850 respectively. h–i, Distribution of enriched mutations for various amino acid substitutions (h) and among different types of base pair substitutions (i) in both saturation synonymous mutations and clinical synonymous mutations. Top score represents the maximum |screen score| value for different epegRNAs corresponding to each mutation in the screen. For c and g, boxplots are depicted as follows: the center line represents the median, the box limits denote the upper and lower quartiles, and whiskers extend to 1.5 times the interquartile range.

Extended Data Fig. 6 Distribution of different key features for synonymous mutations.

a, Distribution of SpliceAI scores for synonymous mutations. b, Distribution of RNA folding energy changes for synonymous mutations. c, Distribution of codon usage frequency changes for synonymous mutations. d, Distribution of conservation scores for synonymous mutations. e, Distribution of gene essentiality (Chronos score) for synonymous mutations in HCT116 cells. f, Distribution of expression levels for synonymous mutations in HCT116 cells.

Extended Data Fig. 7 Validation of synonymous mutations causing aberrant RNA splicing.

a, Schematic depiction of the splicing alterations caused by the RPL11_V132 (GTG > GTC) mutation. The transcript sequence information was obtained by sequencing the cDNA from the experimental and control groups. b, Validation of the effect of the RPL11_V132 (GTG > GTC) mutation on cell proliferation in HCT116 cells. c, Relative mRNA expression levels of RPL11 in the experimental and control groups. The mRNA level of each sample was quantified by real-time qPCR and normalized by GAPDH, and the indicated relative mRNA level was normalized to that of AAVS1-targeting control cells. d, Analysis of editing outcomes for epegRNA targeting RPL11_V132 and controls via genome sequence amplification and NGS. e, Schematic depiction of the splicing alterations caused by the KRAS_G75 (GGG > GGT) mutation. The transcript sequence information was obtained by sequencing the cDNA from the experimental and control groups. f, Validation of the effect of the KRAS_G75 (GGG > GGT) mutation on cell proliferation in HCT116 cells. g, Relative mRNA expression levels of KRAS in the experimental and control groups. The mRNA level of each sample was quantified by real-time qPCR and normalized by GAPDH, and the indicated relative mRNA level was normalized to that of AAVS1-targeting control cells. h, Analysis of editing outcomes for epegRNA targeting KRAS_G75 and controls via genome sequence amplification and NGS. All data are presented as mean ± s.d. (n = 3 biological replicates for cell proliferation assay, n = 3 technical replicates for real-time qPCR). P values were calculated using two-tailed Student’s t test, *P < 0.05, ****P < 0.0001; n.s., not significant.

Extended Data Fig. 8 Editing and the gene expression analysis at the PLK1_S2 site.

a, Analysis of editing outcomes for epegRNA targeting PLK1_S2 and controls via genome sequence amplification and NGS. b, Relative mRNA expression levels of PLK1 in the experimental and control groups. The mRNA level of each sample was quantified by real-time qPCR and normalized by GAPDH, and the indicated relative mRNA level was normalized to that of AAVS1-targeting control cells. Data are presented as mean ± s.d. (n = 3 technical replicates). P values were calculated using two-tailed Student’s t test, n.s., not significant.

Extended Data Fig. 9 Capture and data quality control for DIRECTED-seq.

a, Overview of the process for capturing epegRNA and single-cell transcriptomes. The cDNA of epegRNA was obtained by reverse transcription using primers targeting the constant secondary structure of evopreQ₁, and the cDNA of the entire transcriptome within the cells was obtained by reverse transcription using oligo-dT primers. The figure was adapted from ref. ⁶³. b, Sequencing depth impact epegRNA detection probability and observed gene expression levels. The blue dashed line represents the linear regression line correlating total epegRNAs per cell (y) with total UMIs per cell (x). c–d, Quantile–quantile plots comparing 15 nontargeting epegRNAs (negative control) with 395 epegRNAs across 213 genes (discovery). Genes for the nontargeting tests were randomly selected from the entire gene set. The left-hand panel shows the QQ plot of raw P values, and the right-hand panel shows the QQ plot of –log₁₀ transformed P values. In both panels, the black solid line denotes the theoretical y = x reference under the null hypothesis and the surrounding gray band represents the pointwise 95% confidence interval of that regression.

Extended Data Fig. 10 Comparison of the performance of different prediction models and the screening results in the A549 and KYSE-30 cell lines.

a–b, ROC analysis of DS Finder, SilVA, and CADD on 45 confirmed pathogenic mutations (positive controls) and 1,439 synonymous mutations without phenotypes from our screening (negative controls). The x-axis represents the false positive rate, and the y-axis represents the true positive rate. The red, dark blue, and green curves represent DS Finder, SilVA, and CADD, respectively. The test set in a contains 20 positive controls, excluding the mutations used in the SilVA training set and the test set in b consists of the 25 mutations used in the SilVA training set. c–f, Screening results analysis of three cell lines. Volcano plot illustrating the results of screening for functional synonymous and nonsynonymous mutations affecting cell fitness in A549 (c) and KYSE-30 (e). Blue and red dots denote depleted and enriched epegRNAs, respectively. Analyzing the correlation between three eBARs using zLFC in A549 (d) and KYSE-30 (f), following the same method as in Fig. 1d. The data were filtered according to the log₂ fold change of a specific eBAR and the Pearson correlation of the zLFC of another pair of eBARs was investigated. The y-axis represents the Pearson correlation of the zLFC, while the x-axis indicates the specific threshold that the absolute value of the LFC must surpass. Red, green, and blue denote three individual eBARs: AGCA, CACT, and GCAG, respectively. g, Performance of DS Finder in three cell lines compared with CADD and SilVA. Each point on the graph represents a different dataset, totaling 10.

Supplementary Information

Supplementary Figs. 1–5 and captions for Tables 1–10.

Supplementary Tables 1–10

Supplementary Table 1. The genes involved in the epegRNA library used for this screening. This table provides detailed information on these 3,644 genes and our selection criteria. Supplementary Table 2. Information about the epegRNA library. Sheet 1 provides an overall summary of the library. Sheet 2 contains detailed information, including the target gene, mutation details and sequence information for each epegRNA. Supplementary Table 3. Screening results in HCT116 cells. Supplementary Table 4. Single-cell screening results in HCT116 cells. Supplementary Table 5. Features of synonymous mutations. This table details 75 features used in the synonymous mutation analysis, with comprehensive presentation of the HCT116 screening results. Supplementary Table 6. Clinical prediction results of the DS Finder model and mutation information used for model analysis. Sheet 1 contains prediction scores of colonic disease-associated synonymous mutations by DS Finder. Sheet 2 includes 45 known clinically pathogenic synonymous mutations used for evaluating model prediction performance. Sheet 3 provides same-site controls for these 45 mutations. Sheet 4 contains negative control mutations obtained from screening. Supplementary Table 7. Screening results in A549 cells. Supplementary Table 8. Screening results in KYSE-30 cells. Supplementary Table 9. Sequence information of all epegRNAs used in validation experiments. Supplementary Table 10. All primers used in this study. Each sheet sequentially includes library PCR primers, genomic PCR primers for each validation site, cDNA PCR primers and reverse transcription–qPCR primers.

Cite this article

Niu, X., Tang, W., Liu, Y. et al. Prime editor-based high-throughput screening reveals functional synonymous mutations in human cells.
Nat Biotechnol (2025). https://doi.org/10.1038/s41587-025-02710-z

Download citation

• Received: 10 June 2024

• Accepted: 21 May 2025

• Published: 24 June 2025

• DOI: https://doi.org/10.1038/s41587-025-02710-z