Error-prone PCR is relatively cheap and easy to perform. It uses low-fidelity DNA polymerases to incorporate mistakes during the DNA amplification, and mutation rates can be modified by PCR conditions like different concentrations of manganese chloride and dNTP (Lin-Goerke et al., 1997; Shafikhani et al., 1997). The technique has been widely used in making random mutations for directed evolution experiments (Giver et al., 1998; Moore et al., 2000) and recently for DMS studies (Sarkisyan et al., 2016; Seuma et al., 2021; Faure et al., 2022). However, mutations generated via error-prone PCR are not completely random due to mutation biases of polymerases. For example, Taq polymerase-based mutation rates from A/T are much higher than from C/G (Shafikhani et al., 1997; Wan et al., 1998). Nowadays, error-prone PCR is made easier using commercial kits with mixes of engineered polymerases (Vanhercke et al., 2005), generating reduced biases. However, judging from the DMS data (Faure et al., 2022), mutation biases are only partially removable even with commercial kits. To be noted, error-prone PCR is suitable for generating comprehensive nucleotide-level mutations but not for all possible single amino acid substitutions for each codon. To achieve all possible 19 amino acid substitutions per codon, two consecutive nucleotides of a codon must often be mutated simultaneously. But such a mutation rate will likely hit two or more codons simultaneously, creating a mutation library mixed with single amino acid substitutions and multiple amino acid substitutions.

Another commonly used method is a DMS library with a pool of oligos containing different mutations. Compared to the error-prone PCR, it is more costly but can generate a customized library with fewer biases. Oligonucleotides containing random mutations can be synthesized as a pool of doped oligos (Matteucci and Heyneker, 1983; Araya et al., 2012; Li et al., 2016; Puchta et al., 2016; Li et al., 2019; Wu et al., 2022) or oligos containing NNN triplets (sometimes NNS or NNK) (Hietpas et al., 2012; McLaughlin et al., 2012; Stiffler et al., 2015; Starr et al., 2017; Diss and Lehner, 2018; Hartman et al., 2018; Ahler et al., 2019; Starr et al., 2020; Park et al., 2022), where N represents any of the four nucleotide bases, S for G/C and K for G/T) targeting each codon. This strategy, combined with oligo pool synthesis technology like DropSynth (Plesa et al., 2018), allows construction of user-defined, scalable, and low-cost mutant libraries with comprehensive nucleotide or amino acid substitutions.

These oligos can be designed as doped oligos with each position incorporating a defined percentage of mutations (Starita and Fields, 2015) during the oligo synthesis. The pool of the long mutant oligos (up to 300 nt) can be used as DNA templates. These long oligos need to contain flanking wild-type sequences for primer binding, so they can be amplified and replace the wild-type sequences. On the other hand, short oligos with user-defined mutations or NNN triplets serve as primers. Mutations are introduced to the gene in a manner that is similar to site-directed mutagenesis. The oligos containing NNN triplets are more suited to create mutant libraries covering all possible single amino acid substitutions, while this doped oligo method also targets nucleotide-level mutations as error-prone PCR does. The disadvantage of using oligos with NNN triplets is that it often requires at least two consecutive PCR reactions to generate double amino acid substitutions.

Another popular primer-based method is the nicking mutagenesis (Wrenbeck et al., 2016; Faure et al., 2022), which is developed from a method called Pfunkle (Firnberg and Ostermeier, 2012). Both methods use the circular dsDNA as the template and incorporate mutations using a mix of phosphorylated primers. To remove excessive wild-type template, thymidine in the template is replaced with uracil and degraded after the mutagenesis using the uracil DNA glycosylase and exonuclease III (Exo III) (Firnberg and Ostermeier, 2012). For the same purpose, nicking mutagenesis uses a pair of endonucleases (NtBbvCl and NbBbvCl) that nick one strand of the template dsDNA at a time.

Firstly, a 5’ phosphorylated mutant oligo pool as primers is applied to the NtBbvCI-treated ssDNA template to generate the second-strand DNA with mutations. Then, a second phosphorylated primer without mutations will synthesize the complementary strand for each mutated genetic variant, using the PCR-derived strand with mutations as templates. While other primer-based mutagenesis methods require a pair of primers per mutant, both Pfunkle and nicking mutagenesis requires only one primer per mutant, greatly reducing the cost of oligo synthesis. But to perform Pfunkle or nicking mutagenesis, one needs to ensure high-quality circular DNA and careful design of the primer libraries with a freshly phosphorylated state. Nevertheless, nicking mutagenesis has been rising in popularity for achieving codon-level saturation mutagenesis recently.

For DMS studies aimed at endogenous genetic loci, CRISPR-based technologies (Doudna and Charpentier, 2014; Rees and Liu, 2018) are used. The mutant library can be designed as a sgRNA library targeting intended genetic loci (Wang et al., 2014; Hart et al., 2015; Sadhu et al., 2018; Hanna et al., 2021) or as a donor DNA mutant library for homology-directed repair (HDR) template (Findlay et al., 2014; Sharon et al., 2018; Choudhury et al., 2020; Shen et al., 2022). The donor DNA mutant library can be generated using the methods mentioned above and combined into the backbone flanked by the recombination arms and necessary components. Simultaneous use of two gRNAs also enables multiplexed mutagenesis (Campa et al., 2019). Yet, compared to the ectopic expression of a mutation library, there are much fewer DMS studies performed at the endogenous loci due to additional technical limitations—including sgRNA-dependent uneven editing efficiencies (Wang et al., 2014; Bassalo et al., 2018; Choudhury et al., 2020), low HDR efficiency (Findlay et al., 2014), and high incidences of undetected off-target mutations and editing biases (Cui and Bikard, 2016; Zerbini et al., 2017). The challenges are even more prominent when the mammalian cell genome is the target of saturation mutagenesis (Rees and Liu, 2018).

To sum up, different methods of generating mutation libraries have their own pros and cons (Table 1). The choice of the method should be determined primarily by the purpose. For instance, should the experiment target nucleotide or codon level, single or combinations of mutations? How extensive the mutant library should be, and should the mutations be ectopically expressed or integrated into the genome?

Fitness competition is the most straightforward and economical approach for a high-throughput functional selection. Its logic is that if the gene product is required for cell survival or reproduction, cells carrying functional genetic variants will enrich. In contrast, detrimental variants will deplete over time in a culture medium. As a result, frequency changes of genetic variants can be calculated as fitness scores (Domingo et al., 2018). This strategy does not require special equipment, making it easy to conduct. However, the obtained fitness scores may not necessarily be linearly related to the molecular mechanisms of the mutations, making it complicated to acquire mechanistic insight into the mutational effects on the molecular level (Stein et al., 2019). Besides, marginally detrimental effects on molecular functions may be masked due to the non-linear relationship between fitness and molecular function (Soskine and Tawfik, 2010; Stiffler et al., 2015). It also needs to be noted that mutational effects often alter in different environments (Stiffler et al., 2015; Domingo et al., 2018; Chen et al., 2022).

Therefore, it is essential to select an optimal condition that either reflects the physiological situation best (Starita et al., 2015; Braun et al., 2018; Cantor et al., 2018; Hartman et al., 2018; Staller et al., 2018; Ahler et al., 2019; Matreyek et al., 2020; Mighell et al., 2020) to unveil disease-causing mutations or evolutionary paths of mutations. On the other hand, to easily infer biophysical effects, the fitness assay conditions should be selected to be linearly related to the molecular function (Domingo et al., 2018; Li et al., 2019; Leander et al., 2020; Starr et al., 2020; Faure et al., 2022).

Using the protein stability or binding affinity as a phenotype (Araya et al., 2012; Olson et al., 2014; Starr et al., 2020; Faure et al., 2022) is another widely used method to evaluate mutational effects for a protein-coding gene. This approach can capture essential biophysical effects of mutations and give more mechanistic insights into mutations.

Stability assays often involve tagging the target protein to a reporter, like the green fluorescent protein (GFP) as an indicator of the protein stability (Li et al., 2019; Leander et al., 2020; Matreyek et al., 2020; Park et al., 2022). Cells can be sorted based on the fluorescence levels into several bins, followed by deep sequencing of each sorted subpopulation (Peterman and Levine, 2016; Matreyek et al., 2018). Then, each mutant’s mean fluorescence level is calculated based on the frequencies of each genetic variant in each sorted bin.

In vitro display methods, such as phage display (Araya et al., 2012), yeast display (Klesmith et al., 2017; Starr et al., 2017; Cao et al., 2022), and mRNA display (Olson et al., 2014), detect frequency changes of genetic variants based on the binding affinity of the protein to its ligands. Although the experimental results from such an approach reveal the functional effects of mutations, it does not immediately indicate whether mutations affect the function by changing the protein stability or binding affinity, which is termed biophysical ambiguity hereafter. Nevertheless, it is crucial to resolve the biophysical ambiguity of mutations if we want to predict the combined effects of mutations (Otwinowski et al., 2018; Li and Lehner, 2020) accurately. To overcome this, approaches like combining the binding affinity-based functional assay and the stability-based assay (Starr et al., 2020), or predicting mutants’ biophysical effects by analyzing how mutations combine based on a single assay (Otwinowski et al., 2018) have been shown.

Performing two sets of different experiments (Starr et al., 2020) are often troublesome while predicting folding and binding energy changes based on the protein structures (Capriotti et al., 2005; Schymkowitz et al., 2005; Zhang et al., 2020) is not as accurate as experiment results. Recently, a method called ddPCA (Faure et al., 2022) that uses a relatively simple experimental approach to solve the biophysical ambiguity has been developed, which we will discuss in the following part.

The method called Double Deep Protein-Fragment Complementation Assay (ddPCA) (Faure et al., 2022) is based on the protein-fragment complementation (PCA) assay (Tarassov et al., 2008). In ddPCA, the expression ratios of dihydrofolate reductase (DHFR) fragments are tweaked into two sets so that one assay can detect mutational effects on the stability of the protein (AbundancePCA) while the other detects both stability and protein-protein interactions (BindingPCA) (Figures 1B–F).

BindingPCA uses the traditional PCA method in which two interacting partners are each tagged with interacting partners are each tagged with DHFR[1,2] and DHFR[3] fragments (Figure 1D). Mutations that affect the binding affinity to the ligand and/or the protein stability will reduce the functional DHFR concentration inside the cells and therefore minimize cell survival (fitness) (Figures 1B, C). AbundancePCA, on the other hand, only has one protein-coding gene tagged to one fragment of DHFR (DHFR[3]) while overexpressing the other fragment DHFR[1,2]. This allows the cellular fitness to be solely determined by the limiting concentration of the protein tagged with DHFR fragment (i.e., DHFR [1,2], which reflects the protein stability (Figure 1E). The combination of the BindingPCA and the AbundancePCA serves to determine biophysical effects and resolve biophysical ambiguities (Figure 1F). Compared to other experimental approaches, ddPCA is a much simpler approach to unveil stability and binding affinity of mutations because both AbundancePCA and BindingPCA use the same fitness selection system. ddPCA has been applied to several allosteric proteins and resolves the ‘biophysical ambiguities’, as well as pinpointing allosteric sites systematically (Faure et al., 2022; Weng et al., 2022).

Besides the methods mentioned above, enzyme kinetics can be measured in a dynamic system using microfluidics technology. For instance, the High-Throughput Microfluidic Enzyme Kinetics (HT-MEK) in a DMS experiment allows the systematic investigation of enzymes in an automatically valved microfluidics expression system (Markin et al., 2021).

The most widely used platform in DMS studies has been the Illumina HiSeq platforms for their relatively lower error rates compared to the third-generation sequencing platforms (PacBio or Nanopore sequencing) and the higher cost-effectiveness (cost/base pair) compared to Illumina MiSeq (Shendure et al., 2017; Pfeiffer et al., 2018). However, the HiSeq platforms have a sequencing read length limitation to 300 nt. This makes the identification of long-range epistatic interactions challenging if the mutated region exceeds the length limit. To overcome this, barcoding genetic variants can be applied, first to associate genetic variants with barcodes using Miseq or PacBio sequencing (Puchta et al., 2016; Starr et al., 2020; Wu et al., 2022) and then to perform deep sequencing of the barcodes using HiSeq.

Recently, UMI-based Nanopore sequencing (Zurek et al., 2020) and new circular consensus sequencing (CCS) method using PacBio (Wenger et al., 2019) were developed to increase the accuracies of long-read sequencing platforms to >99.5% (Zurek et al., 2020; Karst et al., 2021). This suggests that Nanopore or PacBio may completely substitute HiSeq for DMS studies in the future.

Experiment design and library preparation for sequencing are utterly important to obtain high-quality data, regardless of the sequencing platforms used. One should always 1) start with a sufficient number of molecules per variant in each mutant library; 2) have multiple independent biological replicates; and 3) avoid experimental bottlenecks, over-sequencing or under-sequencing. Especially, it is essential not to over-sequence as that would impactfully hinder accurate prediction of the variant frequencies (Faure et al., 2020). Read-depths should not be more than that of total expected molecule numbers before generating sequencing libraries but also need to be sufficiently bigger than the expected unique counts of genetic variants.

The phenotype of each genetic variant is often quantified as the normalized relative enrichment scores from the aggregated count data (i.e., after versus before selection) compared to that of the wild type, as shown in Eq. 1 below. E v   = log ⁡ 2 F v , o u t p u t F v , i n p u t − log ⁡ 2 F w t , o u t p u t F w t , i n p u t (1)

F _v,output and F _v,input are the frequencies of a given variant v after selection (output) and before selection (input) respectively, and E _v is the normalized enrichment score of the variant to the wild type. When the phenotype is based on the reporter fluorescence intensity and cell sorting (Starr et al., 2020), functional scores as the mean fluorescence signals are estimated based on the variant counts in each sorted bin and the bin fluorescence parameters (Peterman and Levine, 2016).

To estimate mutants’ phenotypes from the sequencing data accurately and to obtain enough statistical power, correct error detection and propagation are essential. However, it is not a simple task as there are many sources of errors in the typical DMS dataset, including sequencing error, Poisson error, errors from the replicates, and stochastic error (Rubin et al., 2017). Enrich2 (Rubin et al., 2017) and a more recent software DiMSum (Faure et al., 2020) are two good statistical frameworks developed for DMS sequencing data to help users reliably quantify the data and perform error estimation. Enrich2 and DiMSum are both based on the Poisson-based sequence count distribution, and they both integrate the empirical variance into account to estimate the errors. However, the way they handle the empirical variance is different. For instance, Enrich2 takes the mix-effects from the empirical variance, while DiMSum introduced replicate-specific additive and multiplicative modifier terms from empirical variance. Thus, the error estimated from the two models differs (Faure et al., 2020). The only available and direct comparisons between the two statistical software are from Faure and others who developed DiMSum. Based on the 12 datasets examined, Enrich2 and DiMSum performed similarly well on the datasets with little overdispersion, but Enrich2 underestimates errors on the dataset with a lot of overdispersion (Faure et al., 2020). Still, DiMSum is not as widely used as Enrich2 in DMS data analysis, likely because it is still a relatively newly developed pipeline. Either error models from Enrich2 or DiMSum cannot capture systematic errors arising from the experiments, which need to be identified and judged by the researchers using diagnostic plots. After this step, one could select only reliable data based on the error thresholds for further analysis.

Relative mutational effects are often presented in a 2D map, with each x- and y-axis representing mutation position and substitution respectively, and the phenotype in a gradient of filled color as a heatmap (Figure 1). With such a descriptive figure giving an overview of the data, one can quickly judge which positions are more sensitive to mutations and whether certain types of substitutes are more acceptable than others. There are also tools developed to visualize both published and own DMS data. MaveVis, developed as part of the MaveDB (Esposito et al., 2019), allows users to generate heatmaps integrating the protein structural information for each position. It is available for both web-based interfaces and as an R package. Another web-based tool called dms-view (Hilton et al., 2020) can provide a quick exploration of the DMS data to look for specific mutations per site and in the context of protein 3D structure in an interactive manner. Compared to MaveVis, the advantage of dms-view is the integration of the protein 3D structure and logo generation based on mutational effects. However, local users cannot use the software as it is only web-based.

To obtain mechanistic insights into genotype-phenotype maps, machine-learning algorithms (Wang et al., 2014; Hart et al., 2015; Majithia et al., 2016; Klesmith et al., 2017; Rocklin et al., 2017; Weile et al., 2017; Gray et al., 2018; Staller et al., 2018; Song et al., 2021; Faure et al., 2022; Hsu et al., 2022; Leander et al., 2022) and deep learning algorithms (Sarkisyan et al., 2016; Gelman et al., 2021; Faure et al., 2022; Tareen et al., 2022; Vaishnav et al., 2022) are frequently used. Especially, a recently developed python package called MAVE-NN (Tareen et al., 2022) can unveil the one-dimensional latent phenotypes (i.e., a hidden type of biophysical parameter values) that are non-linearly linking genotypes to phenotypes, based on the neural-network algorithm. MAVE-NN has its limitations. For example, it cannot analyse DMS data with only single mutations or mutations affecting more than one type of expected biophysical parameters to reveal the latent phenotypes.