High-throughput leader editing screening identification of functional DNA variation in human genome

Abstract

A gene-leader editing screening platform for identifying functional variations associated with human health and disease, substantially configured to annotate genomes with nucleotide resolutions, accompanied by operable disease prediction and treatment for personalized medical treatment.

Description

[0001] CROSS-REFERENCE TO RELATED APPLICATIONS

[0002] This application claims priority to U.S. Provisional Application No. 63 / 489,418, filed on March 9, 2023, the disclosure of which is hereby incorporated by reference in its entirety for all purposes.

[0003] Government support

[0004] This invention was made with government support from the National Institutes of Health under Grants EY027789, HG009402, AG057497, and DA052713. The government has certain rights in this invention.

[0005] introduction

[0006] Advances in genome sequencing have led to the identification of hundreds of millions of genetic variants in the human population, some of which confer risk for common diseases such as diabetes, neurological disorders, and cancer. 1 A major obstacle to understanding the genetic basis of these complex diseases is the lack of functional annotation of disease-associated variants, particularly because such variants are primarily located in noncoding regions. Increasing evidence suggests that noncoding risk variants may contribute to disease pathogenesis by disrupting gene regulation. Even protein-coding variants identified in individuals with disease are often classified as variants of uncertain significance (VUS). Therefore, to realize the potential of personalized medicine, more precise and higher-throughput functional characterization methods are necessary to elucidate the function of disease-associated variants at base-pair resolution, as well as multiplex analyses across genomic loci.

[0007] Advances in genome editing technologies have enabled the large-scale perturbation and assessment of DNA sequences in desired regions. However, fundamental barriers remain for accurate genome annotation using these methods. For example, CRISPRa, CRISPRi, CRISPR deletion, and CRISPR indel have been used in genetic screening strategies to characterize genes and cis-regulatory regions. 2 Both, but failed to identify the causal variant of the disease. Traditional methods for characterizing DNA variants (SNPs) by homologous recombination knock-in are inefficient and low-throughput. Base editors also have limitations, they introduce specific mutations (C→T, A→G, T→C, G→A), and they have different target efficiencies. 3 Thus, there remain significant gaps in methods for effectively characterizing the role of putative pathogenic variants in human health and disease. Robust, high-throughput methods are urgently needed to accomplish the desired edits at base-pair resolution to better understand the genetic basis of disease.

[0008] Prime editing (PE) is a versatile and precise method of genetic engineering that has been developed to introduce any type of edit, including point mutations, insertions, and deletions 4 . Specifically, PE2 employs a Streptococcus pyogenes Cas9 (SpCas9) H840A nickase and Moloney murine leukemia virus (M-MLV) reverse transcriptase. The spacer in the prime editing guide RNA (pegRNA) guides the Cas9 nickase and M-MLV complex to the target site, while the RT template sequence provides the desired edit information. Thus, both targeting and edit information can be easily programmed in the same pegRNA to perform single nucleotide substitutions, insertions, or deletions. PE3 is a newer iteration of PE that can further improve editing efficiency by using an additional single guide RNA (sgRNA) for nicking cleavage, facilitating replacement of the unedited strand 5 . The ability of prime editors to precisely edit the genome suggests the possibility of high-throughput variant-level genome manipulation. Recently, PE screening was used to identify VUS on the NPC1 locus by transfecting pegRNAs and targeted sequencing of the region based on lysosomal function assays 6 . While transient transfection of PE machinery followed by targeted sequencing of the edited locus is able to identify editing events, its scope is limited to that locus, and therefore, scaling up to assess multiple sites in parallel at a large scale is not feasible. In addition to increasing throughput, improved control of transgene copy number, stable expression of PE machinery, and direct site comparison are desirable. SUMMARY

[0010] The present invention provides methods and compositions for high-throughput prime editing screening to identify functional DNA variants in the human genome.

[0011] In one aspect, the present invention provides a high-throughput prime editing screening to identify functional DNA variants in the human genome configured substantially as disclosed herein.

[0012] In one aspect, the present invention provides a genetic screening platform to identify functional variants associated with human health and disease configured substantially as disclosed herein, and preferably annotates the genome with nucleotide resolution, with actionable disease predictions and therapeutics for personalized medicine.

[0013] In one aspect, the present invention provides a hybrid prime editing screening system configured substantially as disclosed herein, and is capable of characterizing genetic variants at base pair resolution and scale, advancing precise genome annotation for disease risk prediction, diagnosis, and therapeutic target identification.

[0014] In some embodiments, the present application provides a screening, platform, system or method substantially as disclosed herein, wherein the host cells are infected with a lentivirus containing nCas9 and M-MLV reverse transcriptase stable expression cassette (b) to obtain transformed cells with stable expression of nCas9 / M-MLV that enables more efficient pegRNA / ngRNA packaging and lentiviral delivery with higher editing efficiency than co-infection methods (d) to improve PE efficiency and facilitate hybrid screening approaches with lentiviral libraries. Figure 1 b) to obtain transformed cells with stable expression of nCas9 / M-MLV that enables more efficient pegRNA / ngRNA packaging and lentiviral delivery with higher editing efficiency than co-infection methods (d) to improve PE efficiency and facilitate hybrid screening approaches with lentiviral libraries. Figure 1 b) to obtain transformed cells with stable expression of nCas9 / M-MLV that enables more efficient pegRNA / ngRNA packaging and lentiviral delivery with higher editing efficiency than co-infection methods (d) to improve PE efficiency and facilitate hybrid screening approaches with lentiviral libraries.

[0015] In some embodiments, the present application provides a screening, platform, system or method substantially as disclosed herein, wherein the host cells are treated with pegRNA containing a scaffolded RNA motif, for example, these motifs (EvopreQ1, MLV-PK1 and MLV-PK2) at the 3' end of the pegRNA 7-9 b) to obtain transformed cells with stable expression of nCas9 / M-MLV that enables more efficient pegRNA / ngRNA packaging and lentiviral delivery with higher editing efficiency than co-infection methods (d) to improve PE efficiency and facilitate hybrid screening approaches with lentiviral libraries. Figure 5 b) to obtain transformed cells with stable expression of nCas9 / M-MLV that enables more efficient pegRNA / ngRNA packaging and lentiviral delivery with higher editing efficiency than co-infection methods (d) to improve PE efficiency and facilitate hybrid screening approaches with lentiviral libraries. Figure 1 b) to obtain transformed cells with stable expression of nCas9 / M-MLV that enables more efficient pegRNA / ngRNA packaging and lentiviral delivery with higher editing efficiency than co-infection methods (d) to improve PE efficiency and facilitate hybrid screening approaches with lentiviral libraries.

[0016] The present application encompasses all combinations of the specific embodiments described herein as if each combination had been expressly stated. BRIEF DESCRIPTION OF DRAWINGS

[0018] Figure 1a-e. Optimization of PE efficiency in mammalian cells using lentivirus delivery. (a) Different strategies tested for optimization of PE efficiency in MCF7 cell line. Upper: co-infection of three different viruses to deliver PE machinery. Lower: dual pegRNA / ngRNA virus infection of clonal MCF7 line stably expressing nickase Cas9 (nCas9) and Moloney murine leukemia virus reverse transcriptase (M-MLV RT). Also shown are two scaffolds and three different structural RNA motifs tested. (b) Lentiviral constructs used to generate MCF7 clones expressing nCas9 / RT. PuroR, puromycin resistance gene. M-MLV RT, Moloney murine leukemia virus reverse transcriptase. (c) RT-qPCR analysis showing relative expression of nCas9 / RT in different clones, normalized to dCas9 expression of an established CRISPRi iPSC line (yellow). (d) Editing efficiency and indel rate at EMX1 and FANCF loci 2 and 4 weeks after PE installation using two different RNA scaffolds. (e) Improved vectors for expression of pegRNA and ngRNA in PE screens. RTT: reverse transcription template, PBS: primer binding site.

[0019] Figure 2a-h. Functional characterization of the MYC enhancer by PE screening for saturation mutagenesis. (Top) The target enhancer is downstream of MYC. (Bottom) The enhancer region is highly enriched for ATAC-seq, H3K27ac, and H3K4mel ChIP-seq signals. The blue region indicates the region selected for PE screening. (b) (Top) Schematic showing the design of the PE saturation mutagenesis screen at the 716 bp enhancer. Each nucleotide undergoes substitution by three nucleotides with PE. (Middle) Each substitution event is covered by three uniquely designed pegRNA / ngRNA pairs. (Bottom) PE screening workflow. (c) Log2 (fold change) of each substitution per base pair, ordered by its genomic location. Mutations with a significant impact on cell fitness are colored. ATAC-seq signals and conservation scores computed by PhastCons are shown. (d) JARVIS scores for base pairs with different numbers of significant substitutions. Boxplots represent the median, IQR, Q1-1.5xIQR, and Q3+1.5xIQR. Outliers are shown as gray dots. The mean is shown as a red dot. P-values were computed using a two-tailed two-sample t-test. (e) Functional PWMs were created for identifying potential TF binding sites. (f) (Top) ChIP-seq signals for 6 TFs in MCF7. The blue region indicates the core enhancer region. (Bottom) Sequence logo for the core enhancer region generated from the functional PWMs of (e). (g) Matched TF binding sites. (h) (Top) Dense trace showing the nucleotide importance scores for GATA3 and ELF1 binding sites derived by the BPNet model.

[0020] Figure 3a-j. PE screens reveal functional SNPs associated with breast cancer. (a) Overview of Alt and Ref library design. In this design, we included breast cancer-associated variants (SNPs), clinical variants (ClinVar), introduced stop codons (iSTOPs), and non-targeting controls. For each variant, a pegRNA / ngRNA pair was designed to introduce either the Alt or Ref allele. (b) Workflow of PE screens with Alt and Ref libraries. MCF7-nCas9 / RT cells were infected with either lentiviral library. Cells were collected at day 2 and day 32 post-infection. The abundance of pegRNA / ngRNA pairs was deep sequenced in samples collected at day 2 and day 32. The relative effect of each variant was determined based on the relative impact of Alt versus Ref alleles on cell growth. (c) Percentage of significant hits (FDR < 0.05) for Alt / Alt, Het, and Ref / Ref genotypes in MCF7 determined from Alt and Ref PE screens. (d) Functional SNPs (red) with positive or negative effects on cell growth were determined by their relative effects in Alt versus Ref screens. Blue dots represent significant iSTOPs, and black dots represent controls. Red dotted line indicates 0.05 FDR. (e) The absolute effects of identified functional iSTOPs and SNPs were higher than the effects of negative controls (P values calculated by two-tailed two-sample t-test). (f) Genomic distance of SNPs tested at each risk locus relative to the TSS of each gene. Red dots are functional SNPs within the gene body, blue dots are functional SNPs in distal regions, and gray dots are SNPs with no significant effect. (g) Relative enrichment of genomic features of identified functional SNPs (P values calculated by two-tailed Fisher’s exact test). The number of SNPs overlapping with each genomic feature is labeled next to each bar. (h) A Venn diagram showing the number of unique transcription factors (TFs) with differential binding sites centered around functional SNPs. The number of SNPs that alter TF binding sites is also in parentheses. (i, j) Examples of functional SNPs that disrupt TF binding sites. (i) The protective allele of rs12275749 (position shown in f) affects a SMAD3 binding site, and (j) the Alt risk allele of rs66473811 (position shown in f) matches a MAZ binding motif.

[0021] Figure 4a-h. Functional clinical variants identified using PE screening. (a) Clinical functional SNPs (red) with positive or negative effects on cell growth were determined by the relative effect of Alt versus Ref alleles on cell fitness. Blue dots represent significant iSTOPs, and black dots represent negative controls. Red dotted line represents 5% FDR. (b) Effect size of identified functional iSTOPs and clinical variants was greater than that of negative controls (P value calculated by two-tailed two-sample t-test). Boxplot represents median, IQR, Q1 -1.5xIQR, and Q3+1.5xIQR. Red dot represents mean. (c) CADD scores of iSTOPs and clinical variants. (d) Number of identified functional VUSs that caused a shift in each amino acid group. (N, non-polar; P, polar; Pc, positively charged; Nc, negatively charged). (e, f) Lollipop plot of functional VUSs in RAD51C and BARD1, which were mapped to their canonical isoforms. Significant VUSs identified are marked in red. Their effect on cell growth is represented by fold change. (g, h) Lollipop plot of nonsense variants in BRCA1 and BRCA2, which were mapped to their canonical isoforms. Significant hits identified are marked in blue. Their effect on cell growth is represented by fold change.

[0022] Figure 5 a-c. Optimization of PE efficiency in MCF7 cell line.

[0023] (a) Lead editing efficiency and indel rate by co-infection of lentivirus expressing pegRNA, ngRNA, and nCas9 / RT in MCF7 cells. (b) Immunofluorescence staining showing localization of nCas9 / RT (red, FLAG-tagged) in the nucleus (blue, DAPI) in MCF7-nCas9 / RT cells. Scale bar, 1000 μm. (c) Editing efficiency and indel rate by PE using three different structured RNA motifs at the 3’ end of pegRNA 2 and 4 weeks post-infection in MCF7-nCas9 / RT cells.

[0024] Figure 6a-f. Results of enhancer characterization in MCF7 cells and PE screening. (a) CRISPR / Cas9 knockout of the MYC enhancer in MCF7 reduces MYC expression. P-values were calculated using a two-tailed two-sample t-test. (b) Distribution of pegRNA / ngRNA pair read counts in the clonal plasmid library. (c) PCA analysis shows high reproducibility of PE screening between biological replicates. (d) Correlation between PE-induced stop codon positions and their effect sizes. Blue line and P-value were calculated using a generalized additive model. The shaded area represents the 95% confidence interval. (e) (Upper) Locations of the three sensitive base pairs (SBPs) with significant substitutions. (Lower) Cumulative distribution plot of SBPs with three significant substitutions along the MYC enhancer and the formula to calculate the slope of each consecutive bin. (f) Line plot of the slope of each consecutive bin along the MYC enhancer. Red dotted line is the cutoff for significant slopes, which is based on the P-value of 0.05 for slopes derived from Z-scores. Red area is the core enhancer region, derived from the region with a bin slope greater than the cutoff (slope > 0.43).

[0025] Figure 7 a-c. Strategy to prioritize genomic sites and clinical variants in PE screening. (a) Selection of MCF7 growth-related genes from CRISPR / Cas9 knockout and base editing screens in MCF7 cells. (b) Strategy to select breast cancer-related SNPs for lead editing screens. (c) Strategy to select clinical variants for lead editing screens.

[0026] Figure 8 a-f. Quality control and preliminary analysis of PE screening of disease variants. (a) Heatmap with pairwise correlations and hierarchical clustering of read counts for lead editing screens. (b) Pearson correlation between log2(fold change) of iSTOPs in Alt library screens and log2(fold change) of gRNAs in CRISPR / Cas9 knockout screens for each target gene. (c) Volcano plot from Alt library screen results. (d) Volcano plot from Ref library screen results. (e) log2(fold change) of each iSTOP in Alt and Ref library screens. (f) Violin plots showing 5% FDR cutoffs for relative effect analysis to compare Alt and Ref libraries. Numbers above the peaks indicate the relationship of significant data points to total data points in each category using 5% FDR. We used the 5th percentile of P-values from negative controls as an empirical significance threshold to achieve a false discovery rate (FDR) of 5%, which is shown by the red dashed line in d-f.

[0027] Figure 9a-c. Functional VUS examples and their potential consequences. (a) Sequence conservation of RAD51 family proteins. RAD51 family proteins were aligned using MUSCLE. Functional VUS identified in RAD51C by PE screening are labeled. (b) Graph showing the binding region between BARD1 and BRCA1. (c) Alphafold predicted the protein structure of the BARD1 and BRCA1 complex. Two hydrogen bonds were identified between wild-type His36 of BARD1 and Asp96 of BRCA1, but were lost after the BARD1 His36Pro mutation.

[0028] Description of the Invention Embodiments

[0029] Unless otherwise indicated or implied, in these descriptions and throughout the specification, the terms "a" and "an" mean one or more, and the term "or" means and / or. It is to be understood that the examples and embodiments described herein are for illustrative purposes only and that various modifications or changes in light thereof will be apparent to those skilled in the art, and such modifications or changes are intended to fall within the spirit and scope of the application and the appended claims. All publications, patents and patent applications cited herein, including citations therein, are hereby incorporated by reference in their entirety for all purposes.

[0030] Example: High-throughput prime editing screening identifies functional DNA variants in the human genome

[0031] ABSTRACT: Despite great progress in detecting DNA variants associated with human diseases, interpreting their functional impact in a high-throughput and base-pair resolution manner remains challenging. Here, we developed a novel hybrid prime editing screening method that can be applied to characterize thousands of coding and non-coding variants in a single experiment with high reproducibility. To demonstrate its application, we first identified essential nucleotides of a 716 bp MYC enhancer by prime editing-mediated saturation mutagenesis. Next, we applied prime editing screening to functionally characterize 1304 non-coding variants associated with breast cancer and 3699 variants from ClinVar. We found 103 non-coding variants and 156 variants of uncertain significance function through affecting cellular fitness. In summary, we demonstrated a hybrid prime editing screening technique that can characterize genetic variants at base-pair resolution and scale, advancing precise genomic annotation for disease risk prediction, diagnosis, and therapeutic target identification.

[0032] Here, we optimized prime editing (PE) in mammalian cells to enable high-throughput mixed screening of thousands of DNA variants in the human genome via lentiviral delivery. We demonstrate the utility of our novel PE screening method in three different applications, including saturation mutagenesis analysis of a 716 bp enhancer, functional characterization of 1304 breast cancer-associated variants, and assessment of the impact of 3699 clinical variants on cellular fitness. Our results establish the versatility of mixed PE screening in precisely characterizing genetic variants in the human genome.

[0033] Optimization of PE efficiency via lentiviral delivery in mammalian cells

[0034] To enable PE screening with lentiviral delivery, we initially installed PE3 by infecting MCF7 cells with three different viruses: 1) a virus expressing Cas9 (H840A) nickase (nCas9) and Moloney murine leukemia virus reverse transcriptase (M-MLV RT); 2) a virus expressing a pegRNA; 3) a virus expressing a nick sgRNA (ngRNA). Unfortunately, this strategy resulted in PE efficiency below 1% and relatively high indel rates. This is because the efficiency of simultaneously infecting three different viruses in the same cell is low ( Figure 1 a, Figure 5 a)。

[0035] Packaging all PE3 components in the same virus is challenging. To improve PE efficiency and facilitate mixed screening methods with lentiviral libraries, we infected MCF7 cells with a lentivirus containing a stable expression cassette for nCas9 and M-MLV RT ( Figure 1 b). After puromycin selection, we isolated multiple clones and picked the one with the highest nCas9 expression ( Figure 1 c, RT-qPCR, clone #4, Figure 5 b) for subsequent experiments. Stable expression of nCas9 / M-MLV RT enables efficient pegRNA / ngRNA packaging and lentiviral delivery, and the editing efficiency is higher than the co-infection method ( Figure 1 d). To further improve PE efficiency, we evaluated editing efficiency using three different structured RNA motifs (EvopreQ1, MLV-PK1, and MLV-PK2) at the 3’ end of the pegRNA 7-9 Compared to PE using pegRNAs without structured RNA motifs, cells treated with pegRNAs containing the scaffolded RNA motif exhibited consistently higher editing efficiency at both EMX1 and FANCF sites ( Figure 5 c), so we added evopreQ1 in all pegRNA designs for mixed screening. Scaffold 1 5and scaffold 2 10 did not have a significant impact on PE efficiency, suggesting the feasibility of dual delivery of pegRNA and ngRNA from the same viral particle Figure 1 d). All PE experiments in MCF7 cells (MCF7-nCas9 / RT) exhibited relatively low indel rates (0.7% to 1.95%). Therefore, we used MCF7-nCas9 / RT cells and delivered pegRNA with scaffold 1 and ngRNA with scaffold 2 in the same construct using lentivirus Figure 1 e).

[0036] Lead editing enables nucleotide-resolution analysis of enhancer function

[0037] Enhancers can regulate cell type-specific gene expression and are highly enriched for disease-associated variants. Understanding the endogenous function of each nucleotide in an enhancer will reveal key transcription factors that control enhancer activation and aid in the development of better models of gene regulatory networks and in predicting the regulatory effects of disease-associated non-coding variants. To test whether PE screening can quantify the effect of each base in an enhancer, we focused on an MCF7-specific MYC enhancer identified from a CRISPRi screen 11 in MCF7 cells. This enhancer, located 405 kb downstream of MYC, showed enhancer signatures, including open chromatin, H3K27ac, and H3K4mel signals Figure 2 a). Deletion of this enhancer led to 85% downregulation of MYC expression in MCF7 cells, confirming its enhancer activity in regulating MYC expression Figure 6 a). Since MYC downregulation is associated with MCF7 cell survival 12 , we performed a PE-enabled high-throughput saturation mutagenesis screen of this MYC enhancer in MCF7 cells according to the cell survival phenotype Figure 2 b).

[0038] To profile the function of the enhancer at base pair resolution, we designed a library of 6252 pegRNA / ngRNA pairs to generate 2148 single-nucleotide substitutions within the 716 bp MYC enhancer region. Specifically, we changed the original base to three other nucleotides and assessed each event independently three times in the same screen Figure 2b). We also included 94 positive control pegRNA / ngRNA pairs that introduce stop codons (iSTOPs) in the MYC coding region and 400 negative control pegRNA / ngRNA pairs. Of the negative controls, 246 target non-human genomic loci and 154 target the AAVS1 safe harbor site. We then infected MCF7-nCas9 / RT cells with lentiviral libraries expressing these pegRNA / ngRNA pairs Figure 6 b). Two days post-infection, virus-transduced cells were selected with hygromycin for one week and then expanded in regular media for three more weeks. We collected cells at 2 days post-infection and 30 days post-infection, expanded the integrated pegRNA / ngRNA pairs, and determined the relative depletion or enrichment of each pegRNA / ngRNA pair between these two time points by deep sequencing Figure 2 b). We performed 3 such screens Figure 6 c), and used negative controls including pegRNA / ngRNA pairs targeting non-human loci and AAVS1 for data normalization. Using the MAGeCK pipeline 13 Fold change (FC) was calculated for each pegRNA / ngRNA pair between the 2-day and 30-day post-infection samples. As expected, 78% (73 / 94) of the iSTOPs were depleted (log2FC < 0) by 30 days post-infection. The depletion rate of iSTOPs was negatively correlated with their distance to the MYC transcription start site (TSS), consistent with the observation that introducing perturbations at the 5’ end results in higher efficiency of gene knockout 14 Figure 6 d). In addition, two iSTOPs that target the region between the nuclear localization signal (NLS) and the carboxy-terminal domain (CTD) domain were also significantly depleted (amino acid positions 350 and 355) Figure 6 d). The N-terminus of MYC contains its core transcriptional activation domain, which binds multiple partners 15 Those two iSTOPs can create a truncated MYC that is still able to bind cofactors but cannot bind MYC DNA targets, thereby interfering with the function of wild-type MYC and its cofactors.

[0039] ​To investigate the effect of each nucleotide on enhancer function, we defined sensitive base pairs (SBPs) as nucleotides that were affected at least once after substitution to impact cell fitness (FDR < 0.05, |log2FC| > 1). Out of 716 tested base pairs, 334 (46.6%) were SBPs with log2FC < -1, indicating that mutations at these positions decreased enhancer activity and cell fitness. 23.1% (77 / 334) of SBPs were depleted at 30 days after all three substitutions (FDR < 0.05, log2FC < -1). In addition, all tested sequences were not significantly enriched at 30 days, and none increased cell growth phenotype, indicating that perturbation of these sequences completely attenuated enhancer activity. Figure 2 c}.

[0040] Deep learning models have been developed to prioritize non-coding regions and predict their relevance to human diseases. Encouragingly, by JARVIS 16 SBPs with two or more significant substitutions (n = 172) were more deleterious than SBPs with only one significant substitution (n = 162) or non-SBPs (n = 382) Figure 2 d}. This indicates that PE screening was successful in validating the predicted functional sequences by computation. We further established a continuous interval density analysis to detect the change in SBP density along the enhancer to pinpoint SBP enriched regions Figure 6 e and f}. We identified a core enhancer region with high SBP density in the enhancer based on the slope values of the cumulative curve of SBPs with three significant substitutions, as a larger slope value indicates a higher density of SBPs in that region. The core enhancer region was defined with a minimum slope cutoff of 0.43 (P < 0.05 derived from Z-score). The core enhancer region (chr8: 128,142,093-128,142,181, hg38) co-localized with an open chromatin peak. This region contains SBPs with the most significant fold change upon mutation, indicating its strong influence on enhancer activity. Figure 2 c, highlighted in purple). Notably, the core sequence of the enhancer is located next to a highly conserved region Figure 2 c}. This is not surprising as enhancers have experienced rapid evolutionary changes compared to protein-coding sequences 17 .

[0041] Our functional data provided a unique opportunity to compute and construct a position weight matrix (PWM) that captures the functional information of each nucleotide. Using the fold change of each nucleotide substitution from PE screening, we generated a functional PWM Figure 2 e}. We compared our functional PWM with JASPAR, HOCOMOCO, and SwissRegulon databases18-20 Based on the comparison of the selected transcription factor (TF) motifs, 13 TFs with matching motif PWMs were identified Figure 2 g and h). Based on ENCODE ChIP-seq datasets 21 , 5 predicted TFs (GATA3, ELF1, FOXM1, MTA3, and RCOR1) have been shown to bind to MYC enhancers, and Avocado predicted YY1 to bind to this enhancer in MCF7 22 by the ENCODE project Figure 2 f). Furthermore, GATA3 and YY1 are essential cell survival genes in MCF7 23 , which confirmed the usefulness of the saturation mutagenesis approach implemented by our PE to explore enhancer function at base pair resolution. The essential nucleotides for ELF1 and GATA3 binding motifs identified by our PE screen were identical to those estimated by BPNet 24 , further validating the importance of quantifying the role of each nucleotide in the discovery made by our PE screen. Taken together, we demonstrated that hybrid PE screens can be used to elucidate nucleotide resolution functional annotation of non-coding cis-regulatory elements.

[0042] Characterization of breast cancer-associated variants

[0043] Next, we tested the feasibility of characterizing >5000 DNA variants associated with disease at different genetic loci, including non-coding variants from GWAS and variants detected from clinical samples. For GWAS-identified variants, we focused on breast cancer, the most common cancer among women in the United States. To test the feasibility of characterizing DNA variants associated with breast cancer, we used the summary statistics from the largest GWAS to date, which included samples of primarily European ancestry 25 . Candidate genes from the comprehensive fine-mapping effort of this GWAS 26 were selected, which overlapped with genes prioritized by the growth phenotype of the CRISPR screen 23,27 . These genes include CCND1, PSMD6, MYC, UBA52, DYNC1I2, ESR1, MRPS18C, NOL7, EWSR1, BRCA2, and GRHL2 (which were negatively selected in the CRISPR knockout screen) and CUX1, CASP8, and TNFSF10 (which are tumor suppressor genes and were positively selected in the CRISPR knockout screen) (S Figure 7 a). We then selected 1304 single nucleotide polymorphisms (SNPs) within a 500 kbp upstream and downstream range of these genes Figure 7 b), which were previously associated with breast cancer 25related, and it has been suggested that it may act through these genes 26 We also selected 3699 variants from the ClinVar database ( Figure 7 c), 2840 of which were identified from patients undergoing hereditary breast cancer testing 28 To systematically evaluate the effects of variants on cell fitness, we designed two libraries: one introducing a reference allele (Ref library) and the other introducing alternative alleles targeting the selected variants (Alt library) ( Figure 3 a). 250 non-targeting pegRNA / ngRNA pairs were added as negative controls. For the Alt library, 115 pegRNA / ngRNA pairs that introduced stop codons (iSTOP) in 23 MCF7 growth-related genes served as positive controls, while pegRNA / ngRNA pairs that introduced reference sequences were used for these sites in the Ref library. The cloned plasmids were packaged into lentiviral libraries and transduced into MCF7-nCas9 / RT cells. Cells were harvested 2 and 32 days after infection, and the pegRNA / ngRNA pairs were amplified and deeply sequenced ( Figure 3 b). Repeated PE screening using either Ref or Alt libraries (n = 4) was reproducible at the read count level ( Figure 8 a).

[0044] From the Alt library screen, 33.04% (38 / 115) of iSTOPs showed significant cell fitness effects (FDR < 0.05), which is comparable to the 31.8% positive rate of iSTOPs for common essential genes reported in the base editing screen in MCF7 cells. 29 Furthermore, the fold changes of iSTOP were highly correlated with the fold changes of sgRNA targeting the same genes in the MCF7 CRISPR knockout screen. 23 ( Figure 8 b). Compared with day 2, more pegRNA / ngRNA pairs were depleted (FDR < 0.05, Alt PE n = 322 and Ref PE n = 337) than enriched (FDR < 0.05, Alt PE n = 284 = 148 and Ref PE n = 209) in both Alt and Ref PE screening on day 32 (binomial test, P = 4.78x10 -8 , RefPE P = 6.85x10 -16 )( Figure 8c and d). In theory, when designed peg / ngRNA pairs match the wild-type MCF7 genotype, they should have no effect on cell growth. However, it is notable that certain pegRNA pairs that match the wild-type MCF7 genotype exhibit a significantly greater impact on cell growth than expected, and the proportion of significant hits for each genotype group is independent of the initial MCF7 genotype (Chi-squared test for Ref library P = 0.9998, for Alt library P = 0.999, Cochran-Mantel-Haenszel test for Ref and Alt libraries together P = 0.9665). For example, in the Ref library, 11.2% (59 out of 528) of pegRNA for sites with a Ref / Ref MCF7 genotype exhibited a significant depletion, similar to 10.2% (55 out of 540) for heterozygous sites and 7.9% (18 out of 227) for Alt / Alt genotype sites (P = 0.9998, Chi-squared test) Figure 3 c). These changes at sites where alleles are not expected to change suggest that constitutive nCas9 expression has an unintended consequence, similar to CRISPR inhibition (CRISPRi) after the editing machinery is recruited to the target genomic locus 30 To test for potential CRISPRi activity of nCas9 in PE, we compared the results between iSTOPs in the Alt library and the corresponding pegRNA / ngRNA pairs in the Ref library. While pegRNAs in the Ref library exhibited a smaller impact at day 32 than iSTOPs targeting the same loci, they were still depleted at day 32, confirming the unintended consequence due to nCas9 occupying the target genomic loci (P = 0.0001, Wilcoxon rank-sum test) Figure 8 e). Together, we find that prolonged PE expression exhibits an unintended activity similar to CRISPRi, which is a key factor to consider when analyzing lentivirus-mediated PE screens.

[0045] To correct for this unintended PE activity, we corrected each pegRNA / ngRNA pair from the Alt and Ref PE screens for the effect of the SNP allele using DESeq2 31 We compared the FC ratio for each pegRNA / ngRNA pair from the Alt and Ref PE screens. We determined functional SNPs based on the relative impact on cell growth between the Alt and Ref PE. In total, we identified 56 SNPs with the Ref allele and 47 SNPs with the Alt allele that promote cell growth (P < 0.05, empirical significance threshold of 5% to control for Type I error, Figure 8 f, Figure 3 d). As expected, the identified functional SNPs have a smaller effect size than stop codons and a significantly larger effect size than the negative control PE (P = 0.0001, Wilcoxon rank-sum test) Figure 3e). In addition, iSTOPs of genes that promote cell growth, such as MYC and GATA3, were depleted, while iSTOPs of the cell growth inhibitor PTEN were enriched, validating our analytical approach Figure 3 d).

[0046] Because risk variants can be either the Ref or Alt allele, we further annotated functional SNPs based on the gene annotation of breast cancer risk variants. Since most GWAS SNPs can not have causal relationships, we expect that only a fraction of the 1304 tested SNPs will exhibit biological effects. We used CAVIAR to calculate the average likelihood of a variant being causal, and found that the average expectation of a variant being causal is ~8.9% when we assume only one causal variant per linkage disequilibrium (LD) cluster. If we allow for more than one causal variant per LD cluster, the average probability of a variant being causal is ~13.0%. The alternative allele of 50 risk SNPs was promoting cell growth, and the alternative allele of 53 risk SNPs was inhibiting cell growth compared to the reference allele (P < 0.05, two-tailed Fisher's exact test) Figure 3 f). 18.45% (19 / 103) of the functionally validated risk SNPs are located within the risk genes. The remaining SNPs are located in distal regions, with an average distance of 185.8 kb to the TSS of the risk genes Figure 3 f). Except for the BRCA2 locus (where only 2 SNPs were tested), all tested loci contained at least one SNP that had a significant effect on cell growth. Finally, the identified functional SNPs were significantly enriched for active chromatin marks relative to their corresponding genomic background (1 Mbp around the selected cell growth genes), including ATAC-seq, H3K27ac, H3K4mel, and H3K4me3 signals (P < 0.05, two-tailed Fisher's exact test) Figure 3 g).

[0047] To explore the potential mechanisms by which functional SNPs regulate cell adaptive changes, we used the 40 bp region centered at the 103 identified functional SNPs to search for candidate TF binding motifs against the human motif database HOCOMOCO 19 We retrieved 281 and 391 motifs containing the Alt and Ref alleles, respectively (FDR < 0.05 and TF expression > 1 FPKM). After removing redundant motifs for each SNP locus, we identified 90 TF binding sites for 35 unique TFs associated with the cell growth inhibition phenotype (log2FC (Alt / Ref) < 0), and 55 binding sites for 29 unique TFs associated with the pro-cell growth phenotype (log2FC (Alt / Ref) > 0)Figure 3 h). Specifically, the Alt allele (protective allele) rs12275479 (T>C) on the CCND1 locus disrupts a SMAD3 binding motif and is associated with reduced cell growth in PE screening, which is consistent with the TGF -SMAD3 axis reducing MCF7 32 with the number of mammosphere initiating cells in MCF7 cells Figure 3 f and i). In another example, we found that the MAZ binding site of MAZ is affected by the rs66473811 (T>C) Alt allele on the PSMD6 locus. MAZ is a transcription factor that promotes breast cancer cell proliferation by driving tumor-specific expression of the PPARy1 gene and regulating MYC expression 33, 34 which is consistent with the Alt allele as a risk allele Figure 3 f and j). Together, these results support the use of mixed PE screening to functionally characterize GWAS-identified variants.

[0048] Genetic variants detected in clinical samples provide valuable resources for understanding the etiology of human diseases. However, many clinically discovered variants are annotated as variants of uncertain significance (VUS) due to unpredictable functional consequences, even in well-characterized protein-coding genes. To assess the ability of PE screening to functionally annotate VUS using the MCF7 growth phenotype, we designed pegRNA / ngRNA pairs for 2532 VUS, 745 pathogenic variants, and 422 benign variants of 17 genes Figure 7 c). 76.78% of the detected variants were from breast cancer patients. By comparing the relative effect size of each pair of Alt and Ref alleles, we identified 236 functional clinical variants affecting cell growth in 15 genes, including 49 pathogenic variants, 156 VUS, and 31 benign variants Figure 4 a). The average effect size of pathogenic variants, VUS, and benign variants was between the effect size of negative controls and iSTOP Figure 4 b).

[0049] Several computational indices have been used to assess the deleteriousness of variants 35, 36 One such method is CADD, which integrates different genomic annotations into a single quantitative score to assess the relative pathogenicity of human genetic variants 35 iSTOP and pathogenic variants similarly have high CADD scores relative to other categories of variants Figure 4c). VUS and benign variants CADD scores exhibit a broad distribution with median scores much lower than iSTOP and pathogenic variants. Interestingly, in the VUS or benign variant group, the identified functional variants did not have the expected higher CADD scores, suggesting limitations in relying on computational predictions for variant annotation and emphasizing the importance of using functional assays to validate clinical variants, even for variants located in well-studied protein-coding genes. For example, BARD1 (Arg378Ser), a benign variant with a low CADD score (CADD = 4.317), could not be classified as a functional variant. However, according to our PE screening results, this variant exhibited a significant cell growth inhibition effect in MCF7 cells. BARD1 (Arg378Ser) can impair the nuclear localization of the BRCA1 / BARD1 complex and synergize with BARD1 (Pro24Ser) to promote tumor formation in vivo 37 In addition, most of the identified functional VUS are missense variants and about half of the significant VUS we screened changed the amino acid class based on polarity within the same group Figure 4 d), which complicates the determination of their molecular consequences. Our results provide new insights into the potential role of clinical variants in disease pathogenesis through modulating cellular adaptability and provide annotation for previously uncharacterized VUS and benign variants. Functional domains and structural domains are indispensable contributors to protein function. 60% of the identified functional VUS are located within the protein structural domains annotated in the UniProt database 38 Supporting their pathogenicity. For example, we identified 8 VUS in RAD51C Figure 4 e), RAD51C is a cancer susceptibility gene and is also essential for MCF7 survival. Through our PE screening, two variants were associated with reduced cell growth, one variant (Pro21Leu) is located in the RAD51C functional domain (amino acid: 1-126) for Holliday junction processing and the other (Arg366Gln) is located in the NLS region (amino acid: 366-370) Figure 4 e). We also identified functional variants that are not in any annotated domain, including a functional RAD51C VUS (Arg312Gln) associated with a phenotype of reduced MCF7 growth Figure 4 e). Arg312Trp in RAD51C causes homologous recombination deficiency and a phenotype of reduced clonogenicity in MCF10A cells and abolishes RAD51C-RAD51D interaction 39, therefore Arg312Gln may have similar pathogenic consequences on protein function. When comparing the RAD51C sequence with other RAD51 family proteins, we observed that the functional VUS is located at both conserved and non-conserved amino acids ( Figure 9 a), which highlights the challenge of predicting variant function based solely on protein sequence conservation.

[0050] Protein-protein interaction (PPI) is another essential functional activity in many biological processes. In this study, we also identified functional VUSs located in protein binding regions that have the potential to affect PPIs. For example, BARD1 interacts with BRCA1 through its RING domain, and the ubiquitin ligase activity of BRCA1-BARD1 is essential for DNA double-strand break repair. 40, 41 We identified a functional VUS (His36Pro) in the BARD1 RING domain ( Figure 4 f), suggesting that this clinical variant affects the structural consequences of BARD1-BRCA1 heterodimer formation ( Figure 9 b). Consistent with these findings, AlphaFold predicted that the His36Pro variant disrupts the hydrogen bond between His36 in BARD1 and Asp96 in BRCA1 ( Figure 9 c).

[0051] Nonsense mutations generate new stop codons and truncated proteins. Although most nonsense mutations are annotated as pathogenic in ClinVar, the functional consequences of many nonsense mutations remain uncharacterized. 28 In our PE screening, 563 nonsense clinical variants were detected in 13 breast cancer risk genes, of which 38 variants were identified as positive in 7 genes. Notably, 39.47% (15 / 38) showed unexpected phenotypes compared with the cell death knockout phenotypes of these genes. Specifically, a similar number of functional nonsense variants were identified in BRCA1 (n=15) and BRCA2 (n=16) ( Figure 4 g, h); however, 60% (9 / 15) of variants in BRCA1 promoted growth in MCF7 cells, compared with 25% (4 / 16) of variants in BRCA2. After mapping variants in BRCA1 and BRCA2, we noted that all gain-of-function nonsense variants in BRCA1 produced truncated proteins that retained their NLS. These results were confirmed by a different nonsense mutation located at position Q858 downstream of the NLS in BRCA1, which resulted in a truncated BRCA1 with an NLS and increased MCF7 expression. 29 However, for all functional variants identified in BRCA2, their NLS is located at the C-terminus.42 Thus removed from the truncated protein, resulting in loss of BRCA2 nuclear localization. Overall, these results demonstrate the ability of PE screening to functionally characterize certain nonsense mutations.

[0052] Discussion

[0053] We describe a new genome-wide screening method that, by employing and optimizing "search-and-replace" prime editing 5,9 to probe DNA function at base-pair resolution. We demonstrate the success of hybrid prime editing screening in identifying essential nucleotides in the MYC enhancer, functionally characterize 1304 breast cancer-associated risk SNPs, and provide precise annotation of 3699 clinical variants through saturation mutagenic screening. Our study provides a new strategy to elucidate genome function with unprecedented precision and scale. The broad applications demonstrated by this study suggest that hybrid PE screening can significantly enhance the functional characterization toolbox and improve our ability to elucidate the roles of disease-associated variants in the human genome.

[0054] Our analysis indicates that lentiviral PE installation can produce persistent expression of nCas9, pegRNA, and ngRNA, but can lead to unwanted sequence-specific repression similar to CRISPRi. This bias must be corrected to generate precise base-pair resolution annotations. The pegRNA control should be included when assessing the functional impact of a variant, so that other alleles can be introduced at the same locus. Our study normalizes the sequence-specific repression bias by comparing the differential effects of all base-pair substitutions at each locus in the MYC enhancer on cell survival, as well as the differential effects between Alt and Ref alleles in disease variants. Additional improvements can be achieved by controlling the duration of nCas9 expression. For example, doxycycline-inducible nCas9 can be selectively expressed when editing is needed and reversibly turned off afterward. In addition to establishing and optimizing PE screening, we also defined sensitive base pairs (SBPs) and core sequences of MYC enhancer function. We generated a functional PWM for this enhancer by leveraging the effect size of every possible substitution at each base in PE screening. The functional PWM enabled us to precisely predict TF binding sites within the enhancer, providing key annotations to explain MYC activation in MCF7 cells.

[0055] Interpreting the effects of inherited genetic variants will significantly improve our ability to predict disease risk in individuals. However, risk prediction using GWAS data remains limited without substantial functional annotation. In this study, 7.9% of 1304 tested GWAS breast cancer variants, 6.2% of 2532 tested VUS, were identified as significant hits with functions related to MCF7 growth phenotype. Our results demonstrate the feasibility of PE screening to functionally characterize individual variants. It has many applications: for example, PE screening to identify variants associated with differential drug treatment response will help build better predictive models for the unique benefits and risks an individual will derive from a treatment. Using iPSC models, PE screening to read variants directly related to physiological functions, such as lysosomal activity in microglia or synaptic activity in neurons, will reveal functional variants associated with neuropsychiatric diseases. Altogether, our invention provides a functional genomics tool to enable actionable disease prediction, prevention, and treatment necessary for personalized medicine.

[0056] Data availability statement Next-generation sequencing data reported in this study is available from the NCBI Sequence Read Archive database under accession number PRJNA909251. A reviewer link has been created for the bio-project PRJNA909251.

[0057] Methods Cell culture

[0058] MCF7 cells were cultured in Dulbecco's Modified Eagle Medium (DMEM) (Gibco, 10569010) supplemented with 10% fetal bovine serum (FBS) (HyClone, SH30396.03) and passaged with trypsin-EDTA (Gibco, 25200072). All cells were cultured at 37 °C, 5% CO2 and verified to be mycoplasma free using the MycoAlert Mycoplasma Detection Kit (Lonza, LT07-218). Wild-type MCF7 cells were a gift from the Howard Y. Chang lab. The MCF7-nCas9 / RT cell line was generated by lentiviral transduction of cells with an expression cassette expressing a nickase Cas9 (nCas9) Moloney murine leukemia virus reverse transcriptase (M-MLV RT) fusion protein. The infected MCF7 cell pool was treated with puromycin (2.5 µg / ml) for two weeks. Single cells were then sorted by fluorescence-activated cell sorting (FACS) into 96-well plates, one cell per well, to generate clonal MCF7-nCas9 / RT cell lines. The nCas9 / RT expression level in each clone was quantitatively measured by RT-qPCR and normalized to the dCas9 expression level in the WTC11 doxycycline-induced dCas9-KRAB iPSC line 43, 44 .

[0059] Functional characterization of MYC enhancer deletion by CRISPR

[0060] Two sgRNAs were designed to knock out the MCF7 enhancer (chr8: 128,141,747-128,142,627, hg38) (sgl: GAAGTTGTAAGTATAGCGAG, sg2: AGTGCCTGGCACAAGGCAGA). sgRNAs were synthesized in vitro using the Precision gRNA Synthesis Kit (Invitrogen, A29377) according to the manufacturer’s protocol and their concentration was quantified with a Nanodrop. To deliver the genome editing machinery, 100 pmol of Cas9-NLS protein (University of California, Berkeley QB3 MacroLab) and 120 pmo of in vitro synthesized gRNAs were electroporated into 250,000 MCF7 cells using the DN-100 Lonza 4D-Nucleofector program with P3 primary nucleic transfection solution (Lonza, V4XP-3024). Cells were then seeded into 6-well plates and cultured for 2 days before being seeded into 96-well plates to pick single clones. Successfully knocked out clones were identified by genomic PCR with primers forward: CACCAGGACTTGAAGGCAGC, reverse: CACTTCCCAACCTCAGTTTCC. RT-qPCR was used to quantify MYC expression (MYC forward primer: GTCCTCGGATTCTCTGCTCT, reverse primer ATCTTCTTGTTCCTCCTCAGAGTC) and normalized to GAPDH expression levels (GAPDH forward primer: ATTCCATGGCACCGTCAAGG, reverse primer TTCTCCATGGTGGTGAAGACG).

[0061] Cloning of the prime editing plasmid

[0062] To construct the lentiV2-EFla-nCas9 / RT plasmid, we first excised the U6-sgRNA expression cassette from the lentiCRISPR v2 plasmid (Addgene, 52961) by Kpnl and EcoRI double digestion followed by blunt-end ligation. We further replaced the Cas9 expression cassette with the nCas9 / M-MLV-RT expression cassette from the pCMV-PE2 plasmid (Addgene, 132775). The lentiV2-pegRNA and lentiV2-ngRNA plasmids were constructed by replacing the Cas9 and Puromycin sequences in the lentiCRISPR v2 plasmid (Addgene, 52961) with Hygromycin B and EGFP sequences. The RNA motif and sgRNA scaffold were further integrated by Gibson assembly.

[0063] Testing prime editing efficiency

[0064] To assess the prime editing efficiency of EMX1 and FANCF loci, we cloned paired pegRNA / ngRNA into separate vectors. To perform lentiviral co-infection test, we first infected MCF7 cells with EF1a-nCas9 / RT lentivirus, followed by treatment with puromycin (2.5 pg / ml; Sigma-Aldrich, P8833) for 2 weeks to eliminate uninfected cells. Then, EF1a-nCas9 / RT infected cells were seeded in 24-well plates at 12500 cells per well for pegRNA and ngRNA co-infection. 48 hours post-infection, infected cells were treated with hygromycin B (200 pg / ml; Gibco, 10687010), and cells were collected one week post-infection for editing efficiency evaluation. To perform test in MCF7-nCas9 / RT clonal line, we seeded cells in 24-well plates at a density of 12500 cells per well, followed by lentiviral infection (pegRNA-mCherry and ngRNA-EGFP). Two days post-infection, mCherry and EGFP double positive cells were isolated by FACS and cultured. Then, cultured cells were collected at 2 and 4 weeks post-infection for editing efficiency evaluation. 620 genomic DNA was then extracted from each sample using Wizard Genomic DNA Purification Kit (Promega, A1120). The genomic sites of interest were amplified from purified genomic DNA and sequenced on the Illumina NovaSeq 6000 platform. Briefly, for the first round of PCR (PCR1), DNA primers that amplify the genomic sites of interest were used to prepare sequencing libraries. Then, for the second round of PCR (PCR2), DNA primers containing index adapters were used to add these adapters to the PCR1 amplicons. Finally, for the third round of PCR (PCR3), dual-index primers were used to add Illumina indexes to each PCR2 amplicon. Sequencing reads were aligned to the reference sequence using CRISPResso2 45 . For all prime editing efficiency quantification, 21 bp window centered on 1 bp wild type or edited sequence was used to quantify wild type and edited amplicon frequency. Remaining amplicons were categorized as indels. SNP prioritization We selected 14 MCF7 growth-associated genes that overlap with GWAS-identified breast cancer susceptibility genes 26 . For each gene, we selected SNPs from GWAS results from the Breas Cancer Association Consortium 25 . We identified SNPs that are significantly associated genome-wide with P < 1x10-5 , minor allele frequency < 0.02, odds ratio < 0.9 or > 1.2 (roughly representing the upper and lower quartiles of the odds ratio distribution of SNPs that meet the position, P value, and MAF thresholds), these SNPs are located within ± 500 kb of each transcription start site and are associated with breast cancer. We also used a Latino population 46 GWAS results were used to select the ESR1 locus with GWAS P < 1x10 -5 We used the LD Link R package 47 , 641 linkage disequilibrium (LD) clusters were identified among the selected SNPs, with an LD threshold of R 2 > 0.1. Then, we use CAVIAR 48 The most likely causal variants were prioritized, i.e., variants with the causal posterior probability (> 0.1), the highest posterior probability (≤ 0.1), or the largest odds ratio in each haplotype block. We ran CAVIAR twice for each locus, the first time assuming only one causal variant per LD cluster and the second time allowing more than one causal variant per LD cluster. Clinical variant prioritization We retrieved clinical variants from the ClinVar database (accessed 2021-12-25) and retained all single nucleotide variants (SNVs) for the lead editing screening design ( Figure 7 c). We first selected only SNVs that overlapped with genes associated with breast cancer risk and MCF7 growth. Next, we retained only SNVs classified as benign, pathogenic, and of uncertain significance. Furthermore, for SNVs associated with BARD1, BRCA1, BRCA2, RAD51C, RAD51D, and PTEN, we only retained SNVs with more than three submitters, as these genes have thousands of identified variants. Ultimately, our selection criteria yielded 5,310 SNVs, of which we successfully designed pegRNA / ngRNA pairs for 3,699.

[0065] Design and construction of prime editing libraries

[0066] To perform nucleotide-resolution analysis of MYC enhancer function, we first used the PooledDesign-Saturation mutagenesis tool from PrimeDesign. 49PegRNA / ngRNA pairs targeting the 716 bp enhancer region were designed. We optimized the pegRNA / ngRNA pairs based on the proximity of the ngRNA pegRNA (over 50 bp) and the length of the primer binding site (PBS) (close to 14 nt) and redesigned the sequences containing either BsmBI cleavage sites (GAGACG, CGTCTC) or TTTTT. Next, we used GuideScan2 to evaluate the specificity and efficiency of each pegRNA and ngRNA spacer sequence. The low specificity spacer sequences were redesigned to improve the specificity. Finally, three different pegRNA / ngRNA pairs were designed to target the same base pairs, achieving 93.0% (666 / 716) of the substitutions. Each pegRNA / ngRNA pair shared the same pegRNA and sgRNA spacer sequences, only the substitution allele was different in the pegRNA extension sequence. To design the positive control guide, we used pegIT 50 PegRNA / ngRNA pairs were generated that change a single base pair in the MYC coding region to introduce a stop codon. For pegIT 50 For each position suggested, we picked the best pegRNA / ngRNA pair. Based on previous work 51 The AAVS1 locus was chosen as the negative control region for the targeting pegRNA / ngRNA pairs and PrimeDesign 49 Guides were designed as described above. For the non-targeting pegRNA / ngRNA pairs, the pegRNA and ngRNA spacer sequences and pegRNA extension sequences were chosen from the ENCODE non-targeting sgRNA reference dataset (https: / / www.encodeproject.org / files / ENCFF058BPG / ). All pegRNA / ngRNA with non-G 5’-end were added a guanine nucleotide to improve the transcription efficiency of the U6 promoter. We used the following template to join these component sequences: 5’-CTTGGAGAAAAGCCTTGTTT[ngRNA-spacer]GTTTAGAGACG[5nt-random-sequence]CGTCTCACACC[pegRNA-spacer]GTTTTAGAGCTAGAAATAGCAAGTTAAAATAAGGCTAGTCCGTTATCAACTTGAAAAAGTGGCACCGAGTCGGTGC[pegRNA-extension]CCTAACACCGCGGTTC-3’.

[0067] Library oligos for MYC enhancer screening were synthesized by Twist Bioscience and amplified using NEBNext High-Fidelity 2x PCR Master Mix (NEB, M0541L), forward primer: GTGTTTTGAGACTATAAATATCCCTTGGAGAAAAGCCTTGTTT, reverse primer: CTAGTTGGTTTAACGCGTAACTAGATAGAACCGCGGTGTTAGG. To amplify PegRNA / ngRNA library oligos for enhancer saturation mutagenesis pairing, we employed emulsion PCR (ePCR) to reduce recombination of similar amplicons during PCR. Briefly, 96 ePCR reactions of 20 μΐ were performed using 0.01 fmol of mixed oligos and NEBNext High-Fidelity 2x PCR Master Mix (NEB, M0541S). Each 20 μΐ PCR mixture was combined with 40 μΐ oil-surfactant mixture (mineral oil containing 4.5% Span 80 (v / v), 0.4% Tween 80 (v / v), and 0.05% Triton X-100 (v / v)) 52 mineral oil) and vortexed at maximum speed for 5 minutes, centrifuged briefly, and then placed in a PCR machine for amplification. The thermocycler was set to 98°C for 30 s, followed by 26 cycles of 98°C for 10 s, 60°C for 20 s, 72°C for 30 s, followed by 72°C for 5 minutes, and finally 4°C. The ramp rate for each step was 2°C / s. After PCR, individual reactions were combined and purified using QIAQuick PCR Purification Kit (Qiagen, 28104) following previously established guidelines 53The purified PCR product was then treated with Exonuclease I (NEB, M0568L) and purified with 1x AMPure XP beads (Beckman Coulter, A63881). The isolated ePCR products were then inserted into a BsmBI digested LentiV2-mU6-evopreQl vector via Gibson assembly (NEB, E2621L). The assembled products were electroporated into Endura electrocompetent E. coli cells (Biosearch Technologies, 60242) and each library was grown for approximately 4000 individual bacterial colonies. The resulting plasmid DNA was linearized with a BsmbI digest, gel purified, and ligated to a DNA fragment containing the sgRNA scaffold and human U6 promoter using T4 ligase (NEB, M0202M). The resulting library was electroporated into Endura electrocompetent E. coli cells (Biosearch Technologies 60242) and grown as described above. The final plasmid library was extracted using Qiagen EndoFree Plasmid Mega Kit (Qiagen, 12381).

[0068] For SNP and clinical variant screening Alt library, pegRNA / ngRNA pairs were designed using PrimeDesign49. The sequence of 200 bp upstream and downstream of each variant or iSTOP was used as input for PrimeDesign. We used the following parameters to generate initial pegRNA / ngRNA pairs: number of pegRNAs per edit: 10, length of downstream homology: 10 nt, length of PBS: 13 nt, maximum reverse transcription template (RTT) length: 50 nt, number of ngRNAs per pegRNA: 10, nick distance of ngRNA to pegRNA: 50 bp and 75 bp. Next, a guanine nucleotide was added to the 5’ end of all pegRNA / ngRNA pairs with non-G leading nucleotides to improve the transcription efficiency of U6 promoter. The pegRNA / ngRNA pairs containing BsmBI sites (GAGACG, CGTCTC) or TTTTT sequences in pegRNA spacer, ngRNA spacer, or pegRNA extension were excluded. Further selection of pegRNA / ngRNA pairs was made to maximize specificity, efficiency, and ngRNA to pegRNA distance, while minimizing pegRNA to edit distance when multiple pegRNA / ngRNA pairs were available for the same locus. For non-targeting pegRNA / ngRNA pairs, pegRNA spacer, ngRNA spacer, and pegRNA extension sequences were chosen from ENCODE non-targeting sgRNA reference dataset (https: / / www.encodeproject.org / files / ENCFF058BPG / ). To design Ref library, we used the same pegRNA / ngRNA pairs as Alt library but replaced the alternative allele in pegRNA extension with the reference allele. The final oligos followed the following template structure: 5’-CTTGTGGAAAGGACGAAACACC[ngRNA-spacer]GTTTCGAGACG[6nt-random-sequence]CGTCTCTTGTTT[pegRNA-spacer]gttttagagctagaaatagcaagttaaaataaggctagtccgttatcaacttgaaaaagtggcaccgagtcggtgc[pegRNA-extension]TTGACGCGGTTCTATCTAGTTAC-3’.

[0069] Alt and Ref library oligos were synthesized by Twist Bioscience. Alt and Ref plasmid libraries were cloned separately using a two-step cloning approach. First, oligo pools for each library were amplified using NEBNext High-Fidelity 2x PCR Master Mix (NEB, M0541L) and the following primers: forward primer: TCGATTTCTTGGCTTTATATATCTTGTGGAAAGGACGAAACAC, reverse primer: ATTTCTAGTTGGTTTAACGCGTAACTAGATAGAACCGCGTCAA. PCR products were gel excised and column purified (Promega, A9282) and then inserted into BsmBI digested LentiV2-hU6-evopreQl vector by Gibson assembly (NEB, E2621L). Assembly products were electroporated into Endura electrocompetent E. coli cells (Biosearch Technologies, 60242). Each library was grown for ~25 million colonies and then purified using QIAGEN Plasmid Maxi Kit (QIAGEN, 12163). For the second step, plasmid libraries from the first cloning step were linearized by BsmbI digestion, gel purified, and ligated to a DNA fragment containing the sgRNA scaffold and mouse U6 promoter using T4 ligase (NEB, M0202M). Ligated products were electroporated into Endura electrocompetent E. coli cells (Biosearch Technologies, 60242) and each library was grown for ~40 million colonies. Final plasmid libraries were extracted using Qiagen EndoFree Plasmid Mega Kit (Qiagen, 12381).

[0070] Lentivirus production and titration

[0071] To produce lentivirus libraries, we employed the method we previously described 44Briefly, 5 μg of plasmid library was co-transfected with 3 μg of psPAX (Addgene, 12260) and 1 μg of pMD2.G (Addgene, 12259) packaging plasmid into 8 million HEK293T cells plated in a 10 cm dish supplemented with 36 μl of PolyJet (SignaGen Laboratories, SL100688). The medium was changed 12 hours post transfection and then harvested every 24 hours for a total of three times. The harvested viral media was filtered through a Millex-HV 0.45 μm polyvinylidene difluoride filter (Millipore, SLHV033RS) and then further concentrated by centrifugation using a 100000 NMWL (nominal molecular weight limit) Ultra-15 centrifugal filter device (Amicon, UFC910008).

[0072] Lentiviral titers were determined by transducing 400000 cells with increasing volumes (0, 1, 2, 5, 10, 20 and 40 μl) of concentrated virus and polybrene (6 μg / ml; Millipore, TR-1003-G). Forty-eight hours post transduction, cells were dissociated with trypsin-EDTA (0.25%; Gibco, 25200056) and plated in two independent replicates; one set was treated with hygromycin B (200 μg / ml; Gibco, 10687010) for 4 days and the other set was left untreated. Finally, hygromycin resistant cells and control cells were counted to calculate the proportion of infected cells and viral titers.

[0073] Lead editing screen.

[0074] We performed three MYC enhancer PE screens. We transfected MCF7-dCas9 / RT cells with lentiviral libraries at a multiplicity of infection (MOI) of 0.3, with a coverage of 1000 transduced cells per pegRNA / ngRNA pair. After 48 hours, -100 million cells were harvested as a control, and the remaining cells were treated with hygromycin B (200 pg / ml; Gibco, 10687010) for 7 days. After antibiotic selection, cells were maintained in DMEM supplemented with 10% FBS for 30 days, and 100 million cells were collected from the final cell population. We performed four Alt and Ref library screens. We infected -240 million MCF7-nCas9 / RT cells with lentiviral libraries for each replicate of the Alt and Ref screens at an MOI of 0.5, with a coverage of 2000 infected cells per pegRNA / ngRNA pair. Forty-eight hours after infection, one-third of the infected cells from each cell pool were collected as control samples (day 2). The remaining cells were treated with hygromycin B (200 pg / ml; Gibco 10687010) for 7 days and cultured until 32 days after infection (day 32).

[0075] Illumina sequencing library generation

[0076] Genomic DNA was extracted from each sample by cell lysis and digestion [100 mM tris-HCl (pH 8.5), 5 mM EDTA, 200 mM NaCl, 0.2% SDS, and proteinase K (100 pg / ml)], phenol:chloroform (Thermo Fisher Scientific, 17908) extraction, and isopropanol (Thermo Fisher Scientific, BP2618500) precipitation. For MYC enhancer screening, we applied ePCR during library preparation to amplify paired pegRNA / ngRNA sequences in each sample and reduce recombination between similar sequences. Briefly, 400 ng of DNA per reaction was used with NEBNext High-Fidelity 2x PCR Master Mix (NEB, M0541S) and the following primers for 30 20 mI ePCR: Enh-lib-Forward: TCCCTACACGACGCTCTTCCGATCTNNNNNCCTTGGAGAAAAGCCTTGTTT, Enh-lib-Reverse: GGAGTTCAGACGTGTGCTCTTCCGATCTNNNNNGAACCGCGGTGTTAGG. ePCR was performed as previously described to amplify pegRNA / ngRNA pairs from genomic DNA. The thermocycler was set to 98 °C for 30 s, followed by 25 cycles (98 °C 10 s, 60 °C 20 s, 72 °C 1 min), then 72 °C 5 min, and finally hold at 4 °C. The ramp rate was 2 °C / s for each step. After PCR, individual reactions were combined and purified using Ampure XP beads (Beckman Coulter, B23318) following previously established guidelines 53Purified using QIAQuick PCR purification kit (Qiagen 28104). Then, the purified PCR products were treated with exonuclease I (NEB, M0568L) and purified using 1x AMPure XP beads (Beckman Coulter, A63881). The first round of PCR amplicons were used for the second round of PCR to add Illumina adapters and index sequences. For the second round of PCR, we used NEBNext High-Fidelity 2x PCR Master Mix (NEB, M0541S) for 6 ePCR reactions, each containing 0.023 ng of purified DNA. Similar to the first round, the second round of PCR mixtures were prepared and purified. The thermocycler was set at 98°C for 30 s, followed by 12 cycles of 98°C for 10 s, 60°C for 20 s, 72°C for 1 min, then 72°C for 5 min, and finally hold at 4°C. The ramp rate for each step was 2°C / s. For Alt and Ref screening, we amplified the pegRNA / ngRNA pair sequences from each sample using NEBNext High-Fidelity 2x PCR Master Mix (NEB, M0541L) and the following primers: Alt-Ref-lib-forward: TCCCTACACGACGCTCTTCCGATCTNNNNNCTTGTGGAAAGGACGAAACACC, Alt-Ref-lib-reverse: GGAGTTCAGACGTGTGCTCTTCCGATCTNNNNNCGTAACTAGATAGAACCGCGTCAA. Each sample was subjected to 24 50 μΐ PCR reactions, each containing 600 ng of genomic DNA. The individual reactions for each sample were combined and column purified (Promega, A9282). Then, the purified products were amplified by index PCR to add Illumina TruSeq adapters and sample index sequences, with the following primers: index forward: aatgatacggcgaccaccgagatctacac [8 bp index] acactctttccctacacgacgctcttccgatct, index reverse: caagcagaagacggcatacgagat [8 bp index] gtgactggagttcagacgtgtgctcttccgatct. The final libraries were gel purified and subjected to 150 bp paired-end sequencing on the Illumina NovaSeq 6000 platform.

[0077] Data processing and analysis of pilot editing data

[0078] Sequencing libraries were first trimmed with 5 bp random sequences from read 1 and read 2, and low-quality reads were filtered with the fastp tool before formal alignment. To calculate read counts, each pegRNA / ngRNA pair was included if it met the following criteria: (1) read 1 perfectly matched to the sequence containing the 20-21 nt ngRNA spacer and 5 bp flanking sequence; (2) read 2 perfectly matched to the reverse complement of the sequence containing the full pegRNA extension and 5 bp flanking sequence.

[0079] For MYC enhancer PE screen, MAGeCK (0.5.9) pipeline was used 13 The statistical significance and fold change at the sgRNA level for each pegRNA / ngRNA pair, and the statistical significance and fold change at the gene level for each substitution in the cell population relative to the control were evaluated. The pegRNAs targeting non-target and AAVS1 were used as normalized negative controls. To identify the core enhancer region of the MYC enhancer based on the screening results, we first identified base pairs with 3 significant substitutions (FDR < 0.05), and calculated the slope for each interval (moving step = 1 bp, interval size = 30 bp, x-axis: position of each base pair, y-axis: cumulative amount of SBP with 3 significant substitutions). Figure 6 e). The slope was then converted to P-value derived from Z-score accordingly. The core enhancer region was identified by merging overlapping significant stacks (P-value < 0.05).

[0080] For Alt and Ref library screen, zero-read oligos for any sample were removed before proceeding to subsequent analysis. Oligo counts for all samples were input into DESeq2 (1.38.0) 31 and samples with different sequencing depths were normalized using the median of ratios. Then, the normalized read counts for each oligo were modeled as a negative binomial distribution by DESeq2. After that, we used DESeq2 to examine the fold change of each oligo in the Alt and Ref libraries by comparing the data at day 32 and day 2 (design = ~ replicate + condition). We further estimated the relative effect between the reference and alternative alleles by adding an interaction term (design = ~ replicate + condition + allele + condition:allele). The condition refers to the time point of collection (i.e., day 32 or day 2), and the allele refers to the allele class (i.e., Alt or Ref). Finally, Wald test was performed by DESeq2 to calculate P-value. Then, to minimize false positive hits and achieve an empirical FDR of less than 5%, we selected a P-value cutoff corresponding to the fifth percentile of P-values for non-targeting control oligos.

[0081] Motif matrix comparison analysis

[0082] To identify potential transcription factor (TF) binding sites in the target MYC enhancer, we developed a novel method based on motif comparison 54 to directly compare known TF motifs with our base pair resolution functional data. We first computed the log2 (fold change) of each substitution at each base pair with MAGeCK (0.5.9) 13 The log2 (fold change) of the wild-type allele was set to 0. Then, we converted the log2 (fold change) of each substitution to the corresponding fold change value. By normalizing the fold change of each allele per base pair to the sum of all unique allele fold changes per base pair, we further constructed a position weight matrix. We further divided the enhancer sequence into multiple bins of length 5 and 10 base pairs. We only kept bins with information content (IC) above 3 and “N” content less than 10%. Then, we collected all TF motifs (TPM > 10, GSE175204) that are highly expressed in MCF7 cells from the JASPAR, HOCOMOCO, and SwissRegulon databases. Next, we compared the filtered TF motif matrix with the enhancer matrix using Tomtom (P-value < 0.05) to identify potential TF binding sites on the enhancer. Finally, we only kept positive TF motif matches that are at least 95% overlapping with the necessary base pairs of the input sequence (positions with maximum probability > 0.5).

[0083] Predicting base pair contribution to enhancer activity with BPNet

[0084] We trained a convolutional neural network with BPNet that is consistent with published methods 24 to interpret GATA3, ELF1, FOXM1, MTA3, and RCOR1 ChIP-seq data from the ENCODE project. Briefly, the model input is a 1 kb sequence at each ChIP-seq peak site, and the corresponding ChIP-seq control peak is used as the biastrack for training. Regions of chromosome 2 were used as a tuning set, and chromosomes 5, 6, 7, 10, and 14 were used as a test set. Chromosomes X and Y were excluded. The remaining regions of the other chromosomes were used to train the model with default parameters. Once the model for each TF’s ChIP-seq data was obtained, DeepLIFT was used to compute the contribution of each input sequence base pair to the enhancer activity. Finally, TF-MoDISco contribution scores were used to cluster TF motifs and determine merged TF motifs, which were mapped to the input peak region.

[0085] MCF7 genotype analysis

[0086] Sequence Read Archive (SRA) files (paired-end, two reads per locus) for SRR7707725 and SRR7707726 were retrieved from BioProject PRJNA486532. For each run, sequencing reads were individually aligned to the human reference genome hg38 using bwa-mem v.0.7.17. BAM files were then processed using the Picard tools SortSam, MarkDuplicates, and AddOrReplaceReadGroups. Finally, SNPs and indels were called using local haplotype recombination (HaplotypeCaller) using GATK v.4.2.5.0. Single-sample GVCFs (GenotypeGVCFs) from HaplotypeCaller were then jointly genotyped. Finally, genotype concordance between the two runs was verified using CalcMatch v.1.1.2.

[0087] Motif scanning and TF identification of alleles harboring functional breast cancer SNPs

[0088] The 20 bp upstream and downstream sequences of each SNP (Alt and Ref alleles) were used as input sequences for TF motif analysis using FIMO software (version 5.5.0). 55 Comparison with the human TF motif database HOCOMOCO (v11 FULL) 19 To identify matching motifs centered around the SNP region, we used FIMO motif scanning. All FIMO motif scans were performed using default settings. Finally, we selected TFs with binding motifs overlapping the target SNP site (FPKM > 1) (FDR < 0.05, P value < 0.0001).

[0089] Predicting protein structure with AlphaFold

[0090] To explore the effect of the BARD1 His36Pro mutation on the structure of the BARD1 / BRCA1 complex, we used AlphaFold to predict the structures of wild-type BRAD1 / BRCA1 and BARD1 (His36Pro) / BRCA1 complexes. 40 The same amino acid chain used in the determined BARD1 / BRCA1 complex structure (BARD1, residues 26-122; BRCA1, residues 1-103) was used as input for complex structure prediction. The amino acid chains of BARD1 and BRCA1 were imported into the Google Colab version of AlphaFold V2.2.4. 56, 57, which is supported by Python 3 Google Compute Engine. AlphaFold applies a multimer model to respond to the double sequence fill, then searches the genetic database to determine the most suitable multiple sequence alignment (MSA) of the imported sequence, and starts the structure prediction. To avoid stereochemistry violations, all structures are relaxed with the AMBER model (a helper model with energy refinement) and accelerated with GPU. The generated PDB file is imported into UCSF Chimera X 58, 59 for structure visualization. Different colors are given to protein chains to distinguish single chains, and selected amino acid atom structures and hydrogen bonds are drawn out for interaction analysis. Finally, the snapshot function in Chimera X is used to export the real-time rendering of the complex structure at the best visualization angle.

[0091] References

[0092] 1. Taliun, D. et al. Sequencing of 53,831 diverse genomes from the NHLBI TOPMed Program. Nature 590, 290-299 (2021).

[0093] 2. Shalem, O., Sanjana, N.E. & Zhang, F. High-throughput functional genomics using CRISPR-Cas9. Nat Rev Genet 16, 299-311 (2015).

[0094] 3. Anzalone, A.V., Koblan, L.W. & Liu, D.R. Genome editing with CRISPR-Cas nucleases, base editors, transposases and prime editors. Nat Biotechnol 38, 824-844 (2020).

[0095] 4. Chen, P.J. & Liu, D.R. Prime editing for precise and highly versatile genome manipulation. Nat Rev Genet (2022).

[0096] 5. Anzalone, A.V. et al. Search-and-replace genome editing without double-strand breaks or donor DNA. Nature 576, 149-157 (2019).

[0097] 6. Erwood, S. et al. Saturation variant interpretation using CRISPR prime editing. Nat Biotechnol 40, 885-895 (2022).

[0098] 7. Anzalone, A.V., Lin, A.J., Zairis, S., Rabadan, R. & Cornish, V.W. Reprogramming eukaryotic translation with ligand-responsive synthetic RNA switches. Nat Methods 13, 453-458 (2016). 8. Houck-Loomis, B. et al. An equilibrium-dependent retroviral mRNA switch regulates translational recoding. Nature 480, 561-564 (2011).

[0099] 9. Nelson, J.W. et al. Engineered pegRNAs improve prime editing efficiency. Nat Biotechnol 40, 402-410 (2022).

[0100] 10. Dang, Y. et al. Optimizing sgRNA structure to improve CRISPR-Cas9 knockout efficiency. Genome Biol 16, 280 (2015).

[0101] 11. Chen, P.B. et al. Systematic discovery and functional dissection of enhancers needed for cancer cell fitness and proliferation. Cell Rep 41, 111630 (2022).

[0102] 12. Cho, S.W. et al. Promoter of lncRNA Gene PVT1 Is a Tumor-Suppressor DNA Boundary Element. Cell 173, 1398-1412 e1322 (2018).

[0103] 13. Li, W. et al. MAGeCK enables robust identification of essential genes from genome-scale CRISPR / Cas9 knockout screens. Genome Biol 15, 554 (2014).

[0104] 14. Shalem, O. et al. Genome-scale CRISPR-Cas9 knockout screening in human cells. Science 343, 84-87 (2014).

[0105] 15. Baluapuri, A., Wolf, E. & Eilers, M. Target gene-independent functions of MYC oncoproteins. Nat Rev Mol Cell Biol 21, 255-267 (2020).

[0106] 16. Vitsios, D., Dhindsa, R.S., Middleton, L., Gussow, A.B. & Petrovski, S. Prioritizing non-coding regions based on human genomic constraint and sequence context with deep learning. Nat Commun 12, 1504 (2021).

[0107] 17. Villar, D. et al. Enhancer evolution across 20 mammalian species. Cell 160, 554-566 (2015).

[0108] 18. Fornes, O. et al. JASPAR 2020: update of the open-access database of transcription factor binding profiles. Nucleic Acids Res 48, D87-D92 (2020).

[0109] 19. Kulakovskiy, I.V. et al. HOCOMOCO: towards a complete collection of transcription factor binding models for human and mouse via large-scale ChIP-Seq analysis. Nucleic Acids Res 46, D252-D259 (2018).

[0110] 20. Pachkov, M., Balwierz, P.J., Arnold, P., Ozonov, E. & van Nimwegen, E. SwissRegulon, a database of genome-wide annotations of regulatory sites: recent updates. Nucleic Acids Res 41, D214-220 (2013).

[0111] 21. Consortium, E.P. An integrated encyclopedia of DNA elements in the human genome. Nature 489, 57-74 (2012).

[0112] 22. Schreiber, J., Durham, T., Bilmes, J. & Noble, W.S. Avocado: a multi-scale deep tensor factorization method learns a latent representation of the human epigenome. Genome Biol 21, 81 (2020).

[0113] 23. Behan, F.M. et al. Prioritization of cancer therapeutic targets using CRISPR-Cas9 screens. Nature 568, 511-516 (2019).

[0114] 24. Avsec, Z. et al. Base-resolution models of transcription-factor binding reveal soft motif syntax. Nat Genet 53, 354-366 (2021).

[0115] 25. Michailidou, K. et al. Association analysis identifies 65 new breast cancer risk loci. Nature 51, 92-94 (2017).

[0116] 26. Fachal, L. et al. Fine-mapping of 150 breast cancer risk regions identifies 191 likely target genes. Nat Genet 52, 56-73 (2020).

[0117] 27. Hanna, R.E. et al. Massively parallel assessment of human variants with base editor screens. Cell 184, 1064-1080 e1020 (2021).

[0118] 28. Landrum, M.J. et al. ClinVar: improvements to accessing data. Nucleic Acids Res 48, D835-D844 (2020).

[0119] 29. Cuella-Martin, R. et al. Functional interrogation of DNA damage response variants with base editing screens. Cell 184, 1081-1097 e1019 (2021).

[0120] 30. Qi, L.S. et al. Repurposing CRISPR as an RNA-guided platform for sequence-specific control of gene expression. Cell 152, 1173-1183 (2013).

[0121] 31. Love, M.I., Huber, W. & Anders, S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol 15, 550 (2014).

[0122] 32. Bruna, A. et al. TGFbeta induces the formation of tumour-initiating cells in claudinlow breast cancer. Nat Commun 3, 1055 (2012).

[0123] 33. Bossone, S. A., Asselin, C., Patel, A. J. & Marcu, K. B. MAZ, a zinc finger protein, binds to c-MYC and C2 gene sequences regulating transcriptional initiation and termination. Proc Natl Acad Sci U S A 89, 7452-7456 (1992).

[0124] 34. Wang, X. et al. MAZ drives tumor-specific expression of PPARgamma 1 in breast cancer cells. Breast Cancer Res Treat 111, 103-111 (2008).

[0125] 35. Kircher, M. et al. A general framework for estimating the relative pathogenicity of human genetic variants. Nat Genet 46, 310-315 (2014).

[0126] 36. Pollard, K. S., Hubisz, M. J., Rosenbloom, K. R. & Siepel, A. Detection of nonneutral substitution rates on mammalian phylogenies. Genome Res 20, 110-121 (2010).

[0127] 37. Li, W. et al. A synergetic effect of BARD1 mutations on tumorigenesis. Nat Commun 12, 1243 (2021).

[0128] 38. UniProt, C. UniProt: the universal protein knowledgebase in 2021. Nucleic Acids Res 49, D480-D489 (2021).

[0129] 39. Prakash, R. et al. Homologous recombination-deficient mutation cluster in tumor suppressor RAD51C identified by comprehensive analysis of cancer variants. Proc Natl Acad Sci U S A 119, e2202727119 (2022).

[0130] 40. Brzovic, P.S., Rajagopal, P., Hoyt, D.W., King, M.C. & Klevit, R.E. Structure of a BRCA1-BARD1 heterodimeric RING-RING complex. Nat Struct Biol 8, 833-837 (2001).

[0131] 41. Densham, R.M. et al. Human BRCA1-BARD1 ubiquitin ligase activity counteracts chromatin barriers to DNA resection. Nat Struct Mol Biol 23, 647-655 (2016).

[0132] 42. Spain, B.H., Larson, C.J., Shihabuddin, L.S., Gage, F.H. & Verma, I.M. Truncated BRCA2 is cytoplasmic: implications for cancer-linked mutations. Proc Natl Acad Sci U S A 96, 13920-13925 (1999).

[0133] 43. Mandegar, M.A. et al. CRISPR Interference Efficiently Induces Specific and Reversible Gene Silencing in Human iPSCs. Cell Stem Cell 18, 541-553 (2016).

[0134] 44. Ren, X. et al. Parallel characterization of cis-regulatory elements for multiple genes using CRISPRpath. Sci Adv 7, eabi4360 (2021).

[0135] 45. Clement, K. et al. CRISPResso2 provides accurate and rapid genome editing sequence analysis. Nat Biotechnol 37, 224-226 (2019).

[0136] 46. Fejerman, L. et al. Genome-wide association study of breast cancer in Latinas identifies novel protective variants on 6q25. Nat Commun 5, 5260 (2014).

[0137] 47. Machiela, M.J. & Chanock, S.J. LDlink: a web-based application for exploring population-specific haplotype structure and linking correlated alleles of possible functional variants. Bioinformatics 31, 3555-3557 (2015).

[0138] 48. Hormozdiari, F., Kostem, E., Kang, E.Y., Pasaniuc, B. & Eskin, E. Identifying causal variants at loci with multiple signals of association. Genetics 198, 497-508 (2014).

[0139] 49. Hsu, J.Y. et al. PrimeDesign software for rapid and simplified design of prime editing guide RNAs. Nat Commun 12, 1034 (2021).

[0140] 50. Anderson, M.V., Haldrup, J., Thomsen, E.A., Wolff, J.H. & Mikkelsen, J.G. pegIT - a web-based design tool for prime editing. Nucleic Acids Res 49, W505-W509 (2021).

[0141] 51. Chen, C.H. et al. Improved design and analysis of CRISPR knockout screens. Bioinformatics 34, 4095-4101 (2018).

[0142] 52. Williams, R. et al. Amplification of complex gene libraries by emulsion PCR. Nat Methods 3, 545-550 (2006).

[0143] 53. Verma, V., Gupta, A. & Chaudhary, V.K. Emulsion PCR made easy. Biotechniques 69, 421-426 (2020).

[0144] 54. Gupta, S., Stamatoyannopoulos, J.A., Bailey, T.L. & Noble, W.S. Quantifying similarity between motifs. Genome Biol 8, R24 (2007).

[0145] 55. Grant, C.E., Bailey, T.L. & Noble, W.S. FIMO: scanning for occurrences of a given motif. Bioinformatics 27, 1017-1018 (2011).

[0146] 56. Jumper, J. et al. Highly accurate protein structure prediction with AlphaFold. Nature 596, 583-589 (2021).

[0147] 57. Mirdita, M. et al. ColabFold: making protein folding accessible to all. Nat Methods 19, 679-682 (2022).

[0148] 58. Goddard, T.D. et al. UCSF ChimeraX: Meeting modern challenges in visualization and analysis. Protein Sci 27, 14-25 (2018).

[0149] 59. Pettersen, E.F. et al. UCSF ChimeraX: Structure visualization for researchers, educators, and developers. Protein Sci 30, 70-82 (2021).

Claims

1. A high-throughput screening method comprising identifying functional DNA variants in the human genome using pooled prime editing screening.

2. The method according to claim 1, wherein The variants are associated with human health and disease, and the methods further include annotating the genome with nucleotide resolution, leading to actionable disease prediction or treatment for personalized medicine.

3. The method of claim 1, further comprising characterizing genetic variation at base pair resolution and scale to facilitate accurate genome annotation for disease risk prediction, diagnosis, or identification of therapeutic targets.

4. The method according to claim 1, wherein The screen involved dual-pegRNA / sgRNA viral infection of a clonal MCF7 line stably expressing the nickase Cas9 (nCas9) and Moloney murine leukemia virus reverse transcriptase (M-MLV RT).

5. The method according to claim 1, wherein The screen involved MCF7-nCas9 / RT cells and lentiviral delivery of pegRNA with Scaffold 1 and ngRNA with Scaffold 2 in the same construct, as shown in FIG1e .

6. The method of claim 1, wherein The screening involves transfecting host cells with a lentivirus containing a stable expression cassette for nCas9 and M-MLV reverse transcriptase (M-MLVRT), as shown in FIG1b , to obtain transformed cells stably expressing nCas9 / M-MLV RT, which achieve more efficient pegRNA / ngRNA packaging and lentiviral delivery, as well as higher editing efficiency than the co-infection method, to improve PE efficiency and facilitate a mixed screening approach using a lentiviral library.

7. The method of claim 1, wherein The screen involved host cells transfected with pegRNA containing a scaffolding RNA motif at the 3' end of the pegRNA and demonstrated higher editing efficiency at both the EMX1 and FANCF loci compared to using PE without the structured RNA motif.

8. The method of claim 1, wherein The screen comprises host cells transfected with a pegRNA containing a scaffold structured RNA motif at the 3' end of the pegRNA, and exhibits higher editing efficiency at both the EMX1 and FANCF loci compared to using PE without the structured RNA motif, wherein the RNA motif is selected from EvopreQ1, MLV-PK1 and MLV-PK2.

9. The method of claim 1, wherein The screening includes host cells transfected with: (a) a lentivirus containing a stable expression cassette of nCas9 and M-MLV reverse transcriptase (M-MLV RT), as shown in Figure 1b, to obtain transformed cells that stably express nCas9 / M-MLV RT, which achieve more efficient pegRNA / ngRNA packaging and lentiviral delivery, as well as higher editing efficiency than the co-infection method to improve PE efficiency and promote a mixed screening method using a lentiviral library, and (b) a pegRNA that contains a scaffold structure RNA motif at the 3' end of the pegRNA, which exhibits higher editing efficiency at both the EMX1 and FANCF loci compared to using PE without a structured RNA motif.

10. The method of claim 1, wherein The screening includes host cells transfected with: (a) a lentivirus containing a stable expression cassette of nCas9 and M-MLV reverse transcriptase (M-MLV RT), as shown in Figure 1b, to obtain transformed cells that stably express nCas9 / M-MLV RT, which achieve more efficient pegRNA / ngRNA packaging and lentiviral delivery, as well as higher editing efficiency than the co-infection method, to improve PE efficiency and promote a mixed screening method using a lentiviral library, and (b) a pegRNA that contains a scaffold structure RNA motif at the 3' end of the pegRNA, which exhibits higher editing efficiency at both the EMX1 and FANCF loci compared to using PE without a structured RNA motif, wherein the RNA motif is selected from EvopreQ1, MLV-PK1 and MLV-PK2.

11. The method according to claim 1, wherein The screen involved MCF7-nCas9 / RT cells and lentiviral delivery of pegRNA with Scaffold 1 and ngRNA with Scaffold 2 in the same construct, as shown in FIG1e .

Images