Be aware of the allele-specific bias and compositional effects in multi-template PCR

Sample acquisition

The calibration experiment

We extracted DNA using a PowerSoil DNA isolation kit (QIAGEN, Venlo, Netherlands) following the manufacturer’s protocol. To aid our modelling efforts, we tried to minimise all sources of experimental variance as much as possible. To avoid batch effect, all amplicon libraries were constructed in a single run, effectively limiting the number of libraries to the capacity of a T100 thermocycler (Bio Rad, Hercules, CA, USA) we had at our disposal, i.e., 96. Moreover, being restricted to a single amplification run meant opening the thermocycler’s lid and extracting a subset of strip-tubes thereby introducing unwanted disturbances into the amplification process. Therefore, we decided to limit the number of such extraction events and, by extension, the number of cycles to extract. In view of these limitations, we carried out a preliminary qPCR assay mimicking our standard amplicon library preparation procedure to select a range of cycles within the log-linear amplification phase that yielded enough product for subsequent sequencing. The essay were performed with qPCR-HS SYBR (Evrogen, Moscow, Russia) mix accordingly to the manufacturer’s recommendation with universal 16S V4 rRNA primers (F515—GTGCCAGCMGCCGCGGTAA and_R806—GGACTACVSGGGTATCTAAT) (Caporaso et al., 2011). Accordingly to the qPCR essay cycles 22–26 were selected for the deep amplicon sequencing. Further library preparation procedures followed the Illumina MiSeq Reagent Kit Preparation Guide. A single PCR reaction contained 0.15 µl of high-fidelity Encyclo polymerase (Evrogen, Moscow, Russia), 1.5 µl 10X Encyclo buffer (Evrogen, Moscow, Russia), 0.2 µl of 10 mM dNTP (Evrogen, Moscow, Russia), 5.0 pM of universal 16S V4 rRNA primers, 1.0 µl of extracted DNA (approximately 10⁵ templates) and 11.65 µl of PCR water (Evrogen, Moscow, Russia). To avoid temperature fluctuation associated with the edges of the thermocycler, we decided to only fill a centred eight by 12 grid within the thermocycler with 60 tubes covering cycles 22–26 (12 replicates per cycle). To further minimise biases arising from spatial variations in thermodynamic conditions inside the thermocyler, we randomised the placement of strip-tubes (assuming radial symmetry). The libraries were sequenced on an Illumina MiSeq machine using a MiSeq Reagent Kit v3 (600 cycles). Mean number of reads per sample was 49685 (min: 28104, max: 81,330); after the filtration 19,179 reads per samples were remains (min: 9988, max: 31,171). After the filtering, reads were rarefied to 9988 reads.

Data preprocessing

We used DADA2 (Callahan et al., 2016) to infer exact amplicon sequence variants (ASVs) and remove chimera. As per DADA2 recommendations on error-model inference, ASV reconstruction was carried out in separate PCR-cycle groups. To train an IdTaxa taxonomy classifier (Murali, Bhargava & Wright, 2018), we downloaded and preprocessed the SILVA Ref 132 database (Pruesse et al., 2007): we extracted the amplified region, truncated taxonomic descriptions at the highest ambiguous taxonomic rank and removed sequences left without genus-level annotations. To reduce redundancy in the training dataset without mixing different taxa together, we grouped all remaining records by their complete taxonomic identifiers, removed duplicates within the groups using CD-HIT (Fu et al., 2012) and subsampled them without replacement to allow at most 10 representatives per group. The classifier was trained for 30 iterations. We predicted taxonomy and used SEPP (Janssen et al., 2018) to insert ASVs into the reference GreenGenes 13.8 (99%) phylogenetic tree. Amplicon sequence variants that could not be inserted into the tree or were not identified at the genus level with a confidence level ≥ 80% were discarded. To mitigate artifacts of zero-replacement and reduce uncertainty associated with rare amplicon sequence variants, we discarded ASVs that were not observed over 10 times in at least 50% of the libraries. The remaining data were treated by Bayesian-multiplicative zero-replacement with a Dirichlet prior (Martín-Fernańdez et al., 2015).

Isometric log-ratio transform

As frequencies or proportions, HTS data has a compositional nature and requires a particular way of analysis, compositional data analysis (CoDA). A composition is a set of n dependent positive components under the constant sum constraint. To overcome this constraint, a composition should be transformed into the vector of independent components, where the standard analysis for the unconstrained data in real space can be performed. One of these transformations is an isometric log-ratio transform (ILR), which is based on the bipartition strategy developed by Silverman et al. (2017) and Egozcue et al. (2003). Briefly, given a rooted binary tree of n leaves (DNA templates) and n − 1 internal nodes, we can define a sign-matrix Φ of n − 1 rows and n columns such that ${\displaystyle {\Phi }_{ij}=\left\{\begin{array}{cc}-1,\hfill & \text{template}j\text{belongs to the left subclade of node}i\hfill \\ +1,\hfill & \text{template}j\text{belongs to the right subclade of node}i\hfill \\ 0,\hfill & \text{template}j\text{does not descend from node}i.\hfill \end{array}\right.}$

We can now scale rows in matrix Φ and define contrast matrix Ψ of size (n − 1 × n) such that (1)${\displaystyle {\Psi }_{ij}=\left\{\begin{array}{cc}{\Phi }_{ij}\frac{k_i}{n_{i-}},\hfill & {\Phi }_{ij}<0\hfill \\ {\Phi }_{ij}\frac{k_i}{n_{i+}},\hfill & {\Phi }_{ij}>0\hfill \\ 0,\hfill & {\Phi }_{ij}=0\hfill \end{array}\right.}$

where $n_{i+}={\sum }_{j=1}^n\left[{\Phi }_{ij}>0\right]$, $n_{i-}={\sum }_{j=1}^n\left[{\Phi }_{ij}<0\right]$ and $k_i=\sqrt{\frac{n_{i-}{n}_{i+}}{n_{i-}+n_{i+}}}$. The important property of Ψ matrix is that sums of each its row is equal to zero. Let x be a composition, i.e., vector of n positive components, then its isometric log-ratio transform is (2)${\displaystyle \mathrm{ilr}\left(x\right)=\mathrm{ilr}\left(\mathfrak{C}\left(x\right)\right)=log\left(x\right){\Psi }^T=b}$

where ℭ denotes compositional closure (Pawlowsky-Glahn & Buccianti, 2011), b is a vector of so-called balances. Since ILR is an isometry on the Aitchison simplex, we can convert balances back into relative abundances: (3)${\displaystyle x={\mathrm{ilr}}^{-1}\left(b\right)=\mathfrak{C}\left(exp\left(b\Psi \right)\right).}$

PCR amplification model for one sequence

Let z be the amount of an original DNA sequence, $\theta \in \left(0,1\right]$be a constant amplification efficiency associated with the DNA template, and λ ∈ (0, 1] be a constant amplification efficiency of its PCR product. Then, the PCR reaction is expressed as the following recurrence relation capturing a continuous generalisation of the number of PCR products, c, available at time-step (cycle) t: ${\displaystyle c\left(t\right)=\theta z+c\left(t-1\right)+\lambda c\left(t-1\right)=\theta z+\left(\lambda +1\right)c\left(t-1\right),}$

where three terms denote amplified product from the original DNA, remained product from the previous cycle and the amplified product. It should be highlighted, that z and θ variables are present only as a product and are non-identifiable separately. Practically, it means that the initial concentration of a DNA sequence is not detectable in principle, and only an estimate with respect to amplification efficiency, θz, can be obtained. ${\displaystyle c\left(t\right)-c\left(t-1\right)=\left(\lambda +1\right)c\left(t-1\right)-\left(\lambda +1\right)c\left(t-2\right)\Rightarrow c\left(t\right)=\left(\lambda +2\right)c\left(t-1\right)-\left(\lambda +1\right)c\left(t-2\right).}$

We now have a second-degree linear homogeneous recurrence relation, where the initial conditions are c(1) = θz and $c\left(2\right)=\left(\lambda +1\right)\theta z+\theta z$. The relation has the following characteristic equation: ${\displaystyle x^2-\left(\lambda +2\right)x+\left(\lambda +1\right)=0\Rightarrow \left\{\begin{array}{c}{x}_1=\left(\lambda +1\right)\hfill \\ x_2=1.\hfill \end{array}\right.}$

The general solution to this equation is a linear combination of both roots ${\displaystyle c\left(t\right)=\xi {\left(\lambda +1\right)}^t+\zeta 1^t}$

Then, from initial conditions we have ${\displaystyle \left\{\begin{array}{c}c\left(1\right)=\theta z\hfill \\ c\left(2\right)=\left(\lambda +1\right)\theta z+\theta z\hfill \end{array}\right.\Rightarrow \left\{\begin{array}{c}\xi \left(\lambda +1\right)+\zeta =\theta z\hfill \\ \xi {\left(\lambda +1\right)}^2+\zeta =\left(\lambda +1\right)\theta z+\theta z\hfill \end{array}\right.\Rightarrow \left\{\begin{array}{c}\xi =\frac{\theta z}{\lambda }\hfill \\ \zeta =-\frac{\theta z}{\lambda }.\hfill \end{array}\right.}$

Therefore, under our model the number of PCR products associated with a template at cycle t is given by (4)${\displaystyle c\left(t\right)=\frac{\theta z}{\lambda }{\left(\lambda +1\right)}^t-\frac{\theta z}{\lambda }.}$

PCR model for composition data

We model PCR as a discrete-time process parametrised by amplification efficiencies and initial template concentrations, z = (z₁, …, z_n). We distinguish amplification efficiencies associated with original DNA sequences extracted from an environment, θ = (θ₁, …θ_n), and their amplicons, λ = (λ₁, …λ_n), due to stochastic nature of early PCR cycles and differences in template lengths, primer binding-site composition (and, by extension, primer-template complex stability). We assume constant efficiencies (though there is much theoretical and empirical evidence to the contrary (Paliy & Foy, 2011; Chatterjee, Banerjee & Datta, 2012; Kalle, Kubista & Rensing, 2014)) and implicitly restrict our model to the log-linear amplification phase (thereby ignoring potential inter-template competition for substrates at later stages). Under this model (see derivation for Eq. (1)) the number of PCR products associated with template i at cycle t is given by (5)${\displaystyle c_i\left(t\right)=\frac{{\theta }_i{z}_i}{{\lambda }_i}{\left({\lambda }_i+1\right)}^t-\frac{{\theta }_i{z}_i}{{\lambda }_i}}$

where z_i > 0 and ${\lambda }_i,{\theta }_i\in \left(0,1\right]$. Assuming θ_i = λ_i and including the initial concentrations of DNA sequencies, the Expression 5 becomes equivalent to the traditional monoparametric PCR count approximation: (6)${\displaystyle c_i\left(t\right)=z_i{\left({\lambda }_i+1\right)}^t.}$

Assuming large values of t (i.e, PCR cycles >5), we approximate the Expression 5 as follows: (7)${\displaystyle c_i\left(t\right)=\frac{{\theta }_i{z}_i}{{\lambda }_i}{\left({\lambda }_i+1\right)}^t-\frac{{\theta }_i{z}_i}{{\lambda }_i}=\frac{{\theta }_i{z}_i}{{\lambda }_i}{\left({\lambda }_i+1\right)}^t\left(1-{\left(\frac{1}{{\lambda }_i+1}\right)}^t\right)\approx \frac{{\theta }_i{z}_i}{{\lambda }_i}{\left({\lambda }_i+1\right)}^t.}$

While it is impossible to detect absolute amplicon counts, $c_i\left(t\right)$, the observed values after the PCR amplification have a compositional nature. Consequently, we use the isometric log-ratio transform (ILR) and model compositions as balances defined on internal nodes of the amplicon phylogenetic tree (see Eqs. (1) and (2)).

Let $c_i^{\prime }\left(t\right)$ be an observed value of ith template in a sample after t PCR cycles, then the relation between observed and actual amplicon counts is linear, $c_i^{\prime }\left(t\right)=A{c}_i\left(t\right)$, where A is sample-specific factor and the same for all templates of the sample. This factor crucially influences the comparison of amplicon counts between samples, and normalization procedures bring an additional source of bias. To overcome the problem with sample-specific factors, ILR transformation (Eq. (2)) is applied to observed values: ${\displaystyle \mathrm{ilr}\left(c^{\prime }\left(t\right)\right)=log{c}^{\prime }\left(t\right){\Psi }^T=log\left(Ac\left(t\right)\right){\Psi }^T=log\left(c\left(t\right)\right){\Psi }^T+\left(logA,\dots logA\right){\Psi }^T=\left[\text{using property of}\Psi \right]=log\left(c\left(t\right)\right){\Psi }^T=\mathrm{ilr}\left(c\left(t\right)\right).}$

Thus, ILR values of observed amplicon counts equal to ILR values of actual counts. We then substitute Eq. (7) into ILR transformation (8)${\displaystyle \mathrm{ilr}\left(c\left(t\right)\right)=\left(log\left(\lambda +1\right)t-log\lambda +log\theta z\right)\cdot {\Psi }^T=tlog\left(\lambda +1\right){\Psi }^T-log\lambda \cdot {\Psi }^T+log\theta z\cdot {\Psi }^T,}$

where λ = (λ₁, …λ_n), $\theta z=\left({\theta }_1{z}_1,\dots {\theta }_n{z}_n\right)$.

In the Eq. (8), θz is present only once and in the product with Ψ^T, so that the term logθz⋅Ψ^T forms the vector of size (n − 1), and n elements of θz become unidentifiable. This is not surprising due to the fact, that θz reflect biased initial concentrations, which have the compositional nature too. We introduce a = logθz⋅Ψ^T = ilr(θz), and balances of amplicon counts are treated as follows: (9)${\displaystyle b\left(t\right)=tlog\left(\lambda +1\right){\Psi }^T-log\lambda \cdot {\Psi }^T+a.}$

Eq. (9) demonstrate that the growth of PCA products is linear over t in the space of balances. When all parameters of Eq. (8) are estimated ($\hat{\lambda }$ and $\hat{a}$), one can assess relative amplicon counts at any PCR cycle t as follows: (10)${\displaystyle \hat{c}\left(t\right)=\mathfrak{C}\left(\frac{{\hat{z}}_1}{{\hat{\lambda }}_1}\left({\left({\hat{\lambda }}_1+1\right)}^t-1\right),\dots ,\frac{{\hat{z}}_n}{{\hat{\lambda }}_n}\left({\left({\hat{\lambda }}_n+1\right)}^t-1\right)\right),}$

where $\hat{z}={\mathrm{ilr}}^{-1}\left(\hat{a}\right).$

Bayesian inference

Let $x^{tj}=\left(x_1^{tj},\dots ,x_n^{tj}\right)$ be a vector of n measured amplicon counts for jth replicate at cycle t. Then vector of balances of this observation is b^tj = ilr(x^tj). The whole set of observations in our study is ${\left\{b^{tj}\right\}}_{t=\overline{22,26},j=\overline{1,12}}$. We assume that b^tj follows the multivariate Gaussian distribution and take a part in the following hierarchical Bayesian model: ${\displaystyle b^{tj}\sim \mathfrak{N}\left(b\left(t\right),\Sigma \right),b\left(t\right)=tlog\left(\lambda +1\right)\cdot {\Psi }^T-log\lambda \cdot {\Psi }^T+a,\lambda \sim \mathrm{Beta}\left(4,1\right),a\sim \mathfrak{N}\left(0,2\right),\Sigma =\text{diag}\left({\sigma }_i^2\right){\sigma }_i\sim {\mathfrak{N}}_{+}\left(0,1\right).}$

The model was implemented in Python using the probabilistic modelling package PyMC3 (Salvatier, Wiecki & Fonnesbeck, 2015). Parameters were inferred using the No-U-Turn Sampler (Hoffman & Gelman, 2014), a self-tuning variant of Hamiltonian Monte Carlo, over 5000 tuning samples and 20,000 main samples repeated across 4 MCMC chains. We used leave-one-out cross-validation to optimise parameters of prior distributions. We used the Gelman–Rubin statistic and the effective number of samples to monitor convergence and autocorrelation.

Energy estimations

Free energies of secondary structures of 16S rRNA sequencies were calculated at Mfold server (Zuker, 2003; SantaLucia Jr, 1998) with the salt correction (Peyret, 2000). PCR elongation parameters used: temperature 72 °C, salt 50 mM. From a set of alternative structures, we chose those with maximal values of free energies.

Data and code availability

Raw sequencing data were deposited in the NCBI Sequence Read Archive (SRA): BioProject PRJNA545409. All code and metadata required to reproduce the study are openly available on GitHub (https://github.com/arriam-lab2/pcr_bias_publication). We used package ggplot2 (Wickham, 2011) to generate visualisations.

Shift of the microbial community in the PCA amplification

The dataset of amplicon counts consists of 60 observations: five PCR cycles and 12 repeats for each cycle. We performed the amplicon sequence variants picking, and standard PCoA showed the clear shift of the microbial taxonomic structure with the increase of cycles (Fig. 1A). This dynamic is also distinct at the phylum level of the microbial compositions: a rapid expansion of Bacteroides and a much slower growth of Actinobacteria and Proteobacteria at the expense of Firmicutes (Fig. 1B). While, at the resolution of five cycles, these shifts are not very dramatic visually, the PCR bias between the first and 20-th cycles could be significant.

Comparing observed and predicted community dynamics

High-level community dynamics predicted by the PCR amplification model were visually consistent (albeit smoothed out) with observed variations in experimental data (Fig. 1B). On a finer level, experimental log-ratios between relative abundances of individual amplicon sequence variants measured at cycles 26 and 22 were strongly correlated with inferred log-ratios (Pearson r = 0.87, 95% CI = [0.67, 0.95], p = 0). These observations suggest that the amplification model was able to capture and reproduce experimental PCR-induced variations.

Characterising and mitigating amplification bias

We used our model to predicted large-scale community dynamics from cycle 0 to cycle 35 for all ASVs (amplicon sequence variants). PCR bias leads to a dramatic over-representation of Bacteroides and shrinkage of Firmicutes (Fig. 2A). To better visualize the magnitude of amplification biases introduced by a standard amplicon library preparation workflow, we selected amplicon sequence variants (ASVs) classified at the genus level and estimated log-ratio distributions between their relative abundances inferred at cycles 1 (that is, approximated initial template proportions) and 35 (Fig. 2B). Most relative abundances changed by a factor of 2-8. One particularly unfortunate amplicon sequence variant in the Eubacterium eligens group was affected by a factor of 215. The distribution of log-ratios between two distant cycles makes it tempting to view amplification biases as a linear phenomenon.

While the temporal dynamics at the genus level is monotonic, the behaviour of individual ASVs is not that straightforward (Firmicutes’ composition in Fig. 2A). For example, temporal dynamics of individual relative abundances for the top five rare and abundant variants are described by nonlinear curves (Fig. 3A). More importantly, these curves intersect at some PCR cycles and diverge in others, which means the spurious similarity/difference between ASVs abundancies depending on the cycle number. This observation further stresses, that it is absolutely crucial to account for the compositional nature of HTS data.

Investigating putative sources of amplification bias

As we have mentioned in the introduction, it is believed that differences in GC content or secondary structure of DNA templates could be responsible for the amplification bias (Paliy & Foy, 2011; Kalle, Kubista & Rensing, 2014; Kebschull & Zador, 2015; Fan et al., 2019). We also tested the assumption that sequence similarity is a good predictor of effectivity of amplification. But the Mantel test showed no significant correlation between pairwise sequence similarity (measured as Levenstein edit distance) and difference in log-transformed amplification efficiencies.

To evaluate the remaining potential predictors, we fitted a robust linear model for amplification efficiencies λ parametrised by GC content and the energy of secondary structure E: ${\displaystyle log\left(\lambda \right)\sim GC+E.}$

We observed that GC content was not associated with amplification efficiencies, while the energy of secondary structure demonstrated significant association with efficiencies (p-value = 0.0002).