Optimizing illumina next-generation sequencing library preparation for extremely at-biased genomes

Conventional PCR amplification of an AT-rich locus

In order to identify polymerases and PCR conditions that are efficient and tolerant to high AT content templates, we performed conventional PCR using commercial enzymes. The PCR was targeted at a 1217 bp locus (Pf3D7_01:55900-57116) with a high AT-content (84%) and a 540 bp locus (Pf3D7_11:1294982-1295521) with a relatively balanced base composition (positive control) from the P. falciparum 3D7 genome http://www.plasmodb.org[8]. Pf_1 and Pf_11_ama1 oligo pairs (see additional file 1, Table S1) were used to amplify the high AT-rich and positive control fragments respectively. For each enzyme, two sets of PCRs were performed, one using reagents and conditions provided by the manufacturer (standard conditions) and the other deviating only by inclusion of 60 mM TMAC in the reaction mixture (optimized conditions). TMAC is a DNA-binding reagent that has been shown to increase the melting temperature of AT base pairs, thereby improving thermostability of the AT-rich regions during PCR amplification. Although in the presence of 3M TMAC, the thermostability of AT base pairs is equal to GC base pairs, this concentration is also inhibitory to polymerase activity. It has been shown that low concentrations of TMAC (60 mM) improve AT base pair stability and result in increased specificity and yield of PCR products [7]. As shown in Figure 1, most enzymes tested were able to amplify the relatively balanced base composition fragment in the absence of TMAC additive, although with different efficiencies. The presence of TMAC inhibited the activity of most enzymes except for Kapa HiFi, Kapa2G Robust and Platinum pfx, where the additive had a positive effect in amplifying the base-balanced region. At the extreme AT-rich region, nearly all of the enzymes tested were unable to produce detectable amplification either in the presence or absence of TMAC additive. AccuPrime Taq HiFi and Phusion enzymes produced low but detectable amplification of the test region only in the absence of TMAC. However, Kapa HiFi and Kapa2G Robust were the only enzymes that showed efficient amplification of the AT-rich locus in the presence of TMAC additive. Interestingly, the two Kapa enzymes did not produce detectable amplification of this locus in the absence of TMAC, which demonstrates the effect of TMAC on an AT-rich template. PCR amplification of the two test loci were performed at least twice for each enzyme tested and in each case the same results were obtained. Enzymes and PCR conditions that showed some level of amplification of the AT-rich locus were used to amplify genomic libraries for Illumina sequencing. We analyzed the sequence data to assess suitability of each polymerase and conditions for amplification of AT-rich libraries.

PCR amplification of Illumina genomic libraries

Although Phusion polymerase is the standard enzyme for Illumina sequencing library amplification, it produces libraries with poor representation of loci with extreme base composition particularly the AT-rich regions [4]. In our screening above, many enzymes failed to show efficient amplification of the test locus and were clearly intolerant to AT-rich templates. We evaluated the efficiency of genomic library amplification using polymerases that showed some tolerance to the AT-rich template. We tested Phusion, AccuPrime Taq HiFi, Platinum pfx, Kapa2G Robust and Kapa HiFi. We used these enzymes to amplify PE adapter-ligated genomic libraries made from malaria parasite clinical isolate and the reference strain P. falciparum 3D7. In the presence of TMAC, AccuPrime Taq HiFi and Phusion produced undetectable amplifications that were not processed further. However, alternative amplification conditions produced high amplification yields (see additional file 1, Figure S1), that were sequenced and analyzed. All sequence data sets generated were normalized to average coverage before analysis. Two data sets were generated in parallel: One from a lab-cultured parasite (P. falciparum strain 3D7) genomic DNA; and the other one from P. falciparum clinical isolate containing ~53% host DNA contamination. Sequence data from the clinical isolate were filtered to remove host contamination before analysis resulting in lower coverage depths compared to those from 3D7. To analyze and compare the data generated under optimized conditions against standard procedures, we developed analysis metrics that assess four aspects of data quality: Evenness of coverage depth metrics was used to analyze the distribution of coverage depth for each library across the entire genome and plot results of multiple data sets in a single graph for comparative analysis; GC content analysis was used to compute the average GC content for each library and assess base-composition biases; tolerance metrics was used to assess representation of extreme base-composition loci, focusing on selected regions; and fidelity metrics was used to assess enzyme-dependent amplification errors.

Using these metrics, we show that libraries prepared under standard conditions (using Phusion polymerase) produced less even coverage compared to the optimized conditions in both 3D7 and the clinical isolate. As shown in Figure 2A & 2B, up to 30% of the genome was covered at depths of ≤10× for 3D7 and ≤5× for Clinical isolate in Phusion-amplified libraries. In contrast libraries amplified under optimized conditions had most of the genome covered at high depths and only <10% of the genome covered at depths ≤10× and ≤5× for 3D7 and clinical isolate respectively. Libraries prepared under standard conditions also showed strong bias towards the GC-rich templates as shown by a shift in average GC content from ~19.4% (un-amplified reference genome) towards GC neutral base composition (Figure 3A & 3B). Table 1 shows the average GC content for each amplification condition tested. To illustrate this further, we used artemis [9, 10] to view mapped reads on chromosome 1 and assess coverage depth distribution with respect to base composition as depicted by a GC content plot (see Figure 3C and additional file 1, Figure S2 B). We show that under standard conditions, high GC content regions across chromosome 1 including telomeres were over-amplified whereas regions of high AT content were under-amplified. On the other hand Kapa HiFi amplification was close to that of a PCR-free library in coverage of either extremes of base composition.

We also assessed coverage across chromosome 11 and then focused on specific regions of balanced GC content, high AT content and variable base compositions. A box plot showing variation in depth coverage for both P. falciparum strain 3D7 and the clinical isolate is shown on Figure 4i & ii respectively. In these analyses, we observed relatively even coverage depth across the base-balanced region (panel B of Figure 4i & ii) by all the enzymes and conditions tested (this is indicated by the narrow gap between the 25^th and 75^th percentile of base coverage distribution). However, large variations in coverage depths were observed across the AT-rich and the variable base composition loci (Figure 4i & ii panels C and D respectively). Although most enzymes showed more consistency in coverage depth than Phusion, Kapa HiFi produced the best coverage at the AT-rich regions with evenness almost matching that of PCR-free library. At this AT-rich locus in 3D7 (Figure 4i-C), 75th and 25^th percentiles of base coverage depths by both Kapa HiFi libraries produced ~18× and 13× coverage respectively. On the other hand 75th and 25^th percentile of base coverage depth by both Phusion libraries produced ~13× and 3× coverage respectively. The wider gap in coverage depth distribution by Phusion-amplified library indicates uneven and poor coverage of AT-rich regions. For all conditions, coverage of positions 29092-30361 (VAR gene and introns) was generally very low and inconsistent in the clinical isolate (Figure 4(ii) D) compared to 3D7 (Figure 4(i) D), perhaps due to the antigenic variation phenomenon associated with these genes. Overall, we show that our optimized conditions improve the quality of coverage in both clinical and non-clinical samples.

Isothermal exponential amplification of Illumina genomic libraries using RPA technology

High extension temperatures in conventional PCR procedures have been associated with bias and other negative effects on high AT content libraries. We therefore assessed the suitability of an isothermal exponential amplification method provided by the recombinase polymerase amplification (RPA) technology [11]. RPA uses recombinase enzymes which bind oligonucleotide primers (single-stranded DNA) to form nucleoprotein filaments that stimulate searching for homologous sequences in a duplex-DNA template. The primer-recombinase complex is able to attach to the homologous DNA duplex and facilitate strand exchange, eliminating the need to denature the template as in PCR. Disassembly of the recombinase allows strand-displacing DNA polymerase to begin primer extension from the 3' end. Two opposing oligonucleotide primers initiate an isothermal exponential amplification through automatic cyclic repetition of this process. We used RPA technology (TwistDx) with Illumina PE primers to amplify adapter-ligated libraries before sequencing. We obtained a relatively high yield of RPA-amplified library compared to T7 libraries (see additional file 1, Figure S1). Analysis of sequence data obtained from RPA-amplified libraries showed reduced bias towards GC rich templates compared to standard Illumina libraries for both P. falciparum 3D7 and the clinical isolate (Figure 3A & 3B). Coverage of extreme AT base-composition loci was however inferior to Kapa enzymes (Figure 4i & ii panel C). Furthermore, RPA amplified libraries showed relatively high levels of chimeric and duplicate reads (see additional file 1, Table S2).

Linear amplification of Illumina genomic libraries using T7 RNA polymerase

Exponential amplification of sequencing libraries is thought to be the source of biases and artifacts. To avoid this, we explored a recently described linear amplification for deep sequencing (LADS) protocol [12]. LADS uses T7 RNA polymerase to transcribe DNA fragments bearing a T7-promoter sequence at the 5' end. This is an in-vitro transcription system whereby the T7 RNA polymerase transcribes a DNA fragment multiple times producing several RNA transcripts under isothermal conditions and resulting in a linearly amplified template. The resulting RNA material is then converted to cDNA using a reverse transcriptase. We sequenced the T7 polymerase-amplified libraries and analysed the quality of the data generated using our analysis metrics. The overall yield of the library generated using the LADS protocol was very low compared to those amplified exponentially (see additional file 1, Figure S1). However, coverage of extreme base-composition loci was better than standard Illumina libraries. Libraries produced under optimized PCR conditions were generally superior to the LADS protocol (T7) in coverage of AT-rich loci for both P. falciparum 3D7 and the clinical isolate (Figure 4 I & ii panel C). Surprisingly, an over-representation of fragments with high GC content was observed in T7-amplified library when viewing sequence reads aligned to the reference genome using artemis [10] (data not shown). T7 polymerase-amplified libraries were also associated with the highest level of duplicate and chimeric reads (see additional file 1, Table S2).

Effect of TMAC on Illumina genomic library amplification

TMAC has been shown to improve thermostability of AT-base pairs [7]. However, it inhibits the activity of most polymerases. In our analysis, Phusion and AccuPrime Taq HiFi were inhibited by TMAC and produced no library amplification (see additional file 1, Figure S1). On the other hand Platinum pfx and the Kapa enzymes were tolerant to the concentration of TMAC used in this assay (Figure 1). We analyzed and compared libraries prepared with Kapa enzymes in the presence and absence of TMAC. Our analysis indicates that the overall evenness in coverage across the genome was significantly improved in the presence of TMAC (see additional file 1, Figure S3). We observed less bias towards GC rich sequences and an overall improved coverage on AT-rich regions. All amplification data shown for Kapa HiFi, Kapa2G and Platinum pfx in Figures 2, 3 and 4 were performed in the presence of TMAC. The positive effect of TMAC was observed in both P. falciparum 3D7 and the clinical isolate (containing ~53% host contamination).

Fidelity and coverage ranking

Genomic DNA

P. falciparum 3D7 genomic DNA was a gift from Prof Chris Newbold (University of Oxford). Other P. falciparum genomic DNA (PK0076) was a clinical isolate containing ~53% human DNA contamination. The clinical isolate was obtained from the Malaria Genetics Group's Sequencing Sample Repository at the Wellcome Trust Sanger Institute.

Screening for an AT-rich tolerant enzyme

Various commercial enzymes were used to perform a conventional PCR (50 µl total volume) targeted at a 1217 bp locus (Pf3D7_01:55900-57116) with a high AT-content (84%) and a 540 bp locus (Pf3D7_11:1294982-1295521) with a relatively balanced base composition (positive control). Ten ng of P. falciparum 3D7 genomic DNA was used as template for each PCR. Pf_1 and Pf_11_ama1 oligo pairs (see additional file 1, Table S1) were used at a final concentration of 0.3 µM to amplify the high AT-rich and positive control fragments respectively. For each enzyme, two sets of PCRs were performed, one using reagents and conditions provided by the manufacturer (standard conditions) and the other deviating only by inclusion of 60 mM TMAC in the reaction mixture (optimized conditions). Each set of amplification condition was repeated at least twice to confirm the finding. See additional file 1, Table S4 for enzymes used and their catalogue numbers.

Illumina sequencing libraries

Genomic DNA (0.05 - 1 μg in 75 µl TE buffer) from P. falciparum 3D7 or the clinical isolate was sheared for 3 min using a Covaris S2 (Covaris Inc., Woburn, MA).To obtain a fragment-size distribution of ~100 bp to ~300 bp, the settings were: 20% duty cycle, intensity 5, and 200 cycles per burst. Illumina paired-end sequencing libraries were constructed using the NEBNext DNA sample preparation kit (NEB), following standard Illumina sample-preparation protocol with the following slight modifications. End-repair reactions (100 μl) contained 75 µl of sheared DNA sample, 1× T4 DNA-ligase buffer, 1 mM ATP, 0.4 mM dNTPs, 15 units of T4 DNA polymerase, 50 units of Klenow DNA polymerase, 50 units of T4 polynucleotide kinase, and were incubated at 20°C for 45 min. Resulting blunt phosphorylated DNA fragments (end-repaired) were cleaned using the QIAquick PCR purification Kit following the manufacturer's instructions and the DNA eluted in 32 μl of buffer EB. A single A base was added (A-tailing) to the 3' end of the end-repaired DNA fragments in a 50 μl reaction containing 1× Klenow buffer, 0.2 mM dATP, 150 units of Klenow exo- and incubated for 45 min at 37°C. A-tailed DNA fragments were cleaned using QIAquick MiniElute PCR purification Kits following the manufacturer's instructions and DNA eluted in 10 µl EB buffer. Pre-annealed paired-end (PE) adapter oligonucleotides (Illumina) were ligated to the A-tailed fragments in a 50 µl reaction containing 10 µl of DNA sample, 1× Quick T4 DNA-ligase buffer, 10 µl of PE-adapter mixture, 5 µl of Quick T4 DNA ligase (NEB) and incubated at 20°C for 30 min. The ligation reaction was cleaned twice with 1× Agencourt Ampure XP beads (Beckman Coulter). Aliquots were analyzed by running on an Agilent 2100 Bioanalyzer (Agilent Technologies) to determine the size distribution and to check for adapter contamination. Most amplified or non-amplified libraries were prepared twice (starting from same non-sheared genomic DNA) and each technical replicate sequenced independently before data analysis. All replicate data sets were pooled after normalizing and before analysis with the exception of coverage analysis on specific genome loci (presented as box plots, Figure 4i & 4ii) where each replicate was analyzed and shown independently.

Illumina PCR library amplification

All PCR amplifications were performed with an MJ Research Thermo Cycler PTC-225. Illumina PE 1.0 and 2.0 primers or PE 1.0- and 2.0-derived indexing primers were used to amplify 10 ng adapter-ligated library fragments by PCR. From a single pool of adapter-ligated library, aliquots were amplified using standard Illumina PCR reagents or using various alternative polymerases under optimized PCR conditions. Standard Illumina PCR (50 μl) with Phusion polymerase (Thermo) contained 1× Phusion DNA polymerase master mix and 0.4 μM of each primer pair and was amplified with thermocycling conditions of 1 min at 98°C for the initial denaturation, followed by 12 cycles of [10 s at 98°C, 30 s at 65°C, 30 s at 72°C] and a final extension for 5 min at 72°C. PCR (50 μl) with Kapa HiFi (KAPA Biosystems, South Africa) contained 1× Kapa HiFi buffer (containing TMAC) or 1× Kapa buffer B, 0.3 mM of each dNTP, 0.4 μM of each primer pair, 1 unit of Kapa HiFi and was amplified with thermocycling conditions of 1 min at 98°C for the initial denaturation followed by 12 cycles of [10 s at 98°C, 1 min at 65°C] and a final extension for 5 min at 65°C. PCR (50 μl) with Kapa2G Robust (KAPA Biosystems, South Africa) contained 1× Kapa2G buffer B, 0.3 mM dNTPs mix, 0.4 μM of each primer pair, 0.4 units of Kapa2G Robust and was amplified in the presence or absence of 60 mM TMAC with the same thermocycling conditions as Phusion above. PCR (50 μl) with Platinum pfx (Invitrogen) contained 1× Platinum pfx buffer, 0.3 mM dNTPs mix, 0.4 μM of each primer pair, 1 mM MgCl₂, 0.8 units of Platinum pfx and was amplified in the presence or absence of 60 mM TMAC with thermocycling conditions of 2 min at 94°C for the initial denaturation followed by 12 cycles of [15 s at 94°C, 30 s at 60°C, 30 s at 65°C] and a final extension for 5 min at 65°C. PCR (50 μl) with AccuPrime Taq HiFi (Invitrogen) contained 1× AccuPrime buffer II, 0.4 μM of each primer pair, 1 unit of AccuPrime Taq HiFi and was amplified in the presence or absence of 60 mM TMAC with the thermocycling conditions used for Platinum pfx. PCR products were purified with 1× Agencourt Ampure XP beads (Beckman Coulter) and eluted in EB buffer (Qiagen). PCR products were analysed by running on an Agilent 2100 Bioanalyzer (Agilent Technologies) to assess amplification. Only successful amplifications were processed for sequencing.

Illumina isothermal library amplification using RPA technology

We used recombinase polymerase amplification (RPA) technology as provided in the TwistAmp™ DNA-amplification kit [11] to amplify 10 ng adapter-ligated library fragments. Illumina PE 1.0 and 2.0 primers or PE 1.0- and 2.0-derived indexing primers (each in 10 μM stock) were used for amplification following kit manufacturer's instructions. Amplifications were performed twice to generate technical replicates. For each sample, a rehydration solution was prepared by mixing 2.4 μl of each primer pair, 29.5 μl of rehydration buffer and 13.2 μl of template library solution. After vortexing and brief centrifugation, the 47.5 μl of rehydration solution was transferred to a freeze-dried RPA reaction complex pellet and mixed by pipetting up and down until complete resuspension was achieved. To initiate the reaction, 2.5 μl of 280 mM magnesium acetate was added and mixed well by vortexing followed by brief centrifugation. The reaction was allowed to continue in an incubator block at an isothermal temperature of 37°C for 40 min. Amplification products were purified with 1× Agencourt Ampure XP beads (Beckman Coulter) and eluted in EB buffer (Qiagen).

Illumina linear library amplification using T7 RNA polymerase

Linear amplification for deep sequencing (LADS) [12] was used to amplify end-repaired and A-tailed genomic DNA library fragments. Starting with 200 ng of genomic DNA, ~30 ng of A-tailed library fragments were ligated to Adapter-A and Adapter-B, previously prepared by annealing adapter-A and adapter-B forward and reverse oligos (see additional file 1, Table S1) respectively. Adapter-ligated libraries were cleaned twice with 1× Agencourt Ampure XP to remove non-ligated adapters and adapter dimers. Alternatively, adapter-ligated library was electrophoresed in a 2% agarose gel alongside a DNA ladder and fragment bands of 250 to 600 bp were cut and purified at room temperature using the Qiaquick gel extraction kit (Qiagen). The purified adapter-ligated library (20 ng) was used to perform in vitro transcription at 37°C overnight with reagents from the MEGAscript T7 kit (Ambion). The protocol described in [12] was followed to generate T7-amplified libraries ready for Illumina cluster generation and sequencing. The procedure was repeated twice, each time starting from same genomic DNA, to generate technical replicates

PCR-free libraries

Genomic DNA (1-2 μg) was sheared, end-repaired and A-tailed as described under "Illumina sequencing libraries". Pre-annealed paired-end Illumina No-PCR adapters were ligated using the same ligation procedure described under "Illumina sequencing libraries". Resulting adapter-ligated fragments were loaded alongside a low molecular weight ladder (NEB) on a TBE 2% agarose gel. Fragment bands of 350 to 500 bp were cut and purified using a Qiaquick gel extraction kit. Gel purification was performed at room temperate and DNA eluted in EB buffer.

Library quantification by qPCR

Pre-PCR libraries were quantified by qPCR using a previously quantified and sequenced library in dilutions containing 100 ng, 10 ng and 1 ng as standard. The qPCR reaction (25 μl) contained 12.5 μl of 2× Sybr Green PCR Master Mix (Applied Biosystems), 0.75 μl of 10 μM Pre-PCR qPCR top and bottom primers (see additional file 1, Table S1), 2.5 μl of library sample (pre-diluted 1000 fold) and 8.5 μl of PCR-quality water. The reaction was set in triplicates and run on a StepOne™ Real-Time PCR System (Applied Biosystems) with a thermocycling protocol of 10 min at 95°C for the initial denaturation followed by 40 cycles of [30 s at 95°C and 1 min at 60°C]. Prior to cluster generation, amplified or no-PCR libraries were quantified by qPCR using the SYBR fast Illumina library quantification kit (KAPA Biosystems). The qPCR was done following instructions provided by the kit manufacturer and run on the StepOne™ Real-Time PCR System (Applied Biosystems). The thermocycling protocol was 5 min at 95°C for the initial denaturation followed by 35 cycles of [30 s at 95°C and 45 s at 60°C].

Data processing

Replicate data sets were analyzed independently by mapping sequence reads to the P. falciparum 3D7 reference genome V2.1.5 ftp://ftp.sanger.ac.uk/pub/pathogens/Plasmodium/falciparum/3D7/3D7.version2.1.5/ using BWA [13]. SAMtools [14] was used to generate coverage data from the BWA pileup mapping output. To compare the quality of sequence data from different library preparation methods, we developed four analysis metrics. Evenness of coverage depth metrics - comparing the overall representation and depth across the entire genome; GC content analysis - to detect base composition biases in amplification through comparison with a theoretical GC content derived from shredded reads of the reference genome [6]; Tolerance metrics - to assess representation of extreme base-composition loci, focusing on selected genomic regions; and Fidelity metrics - assessing enzyme-dependent errors.

Genome coverage

The value lci (0.5) that gives a measurement of the coverage below half of the average depth in the distribution was used to compare evenness of coverage for each data set.

We counted the number of bases in the genome with zero coverage (Cov = 0) and less than 5× coverage (Cov<5×). We used SAMtools to generate coverage plots and bash/awk scripts for coverage counting. Tolerance to regions of extreme base composition was analyzed after normalizing data sets by overall average coverage depth. Coverage distribution on selected regions (on chromosome 11) are presented as a boxplot generated in R [15]http://www.r-project.org/.

Fidelity

Most polymerases introduce errors during base incorporation in regions of low complexity. To assess these forced errors (or template-dependent errors), we compared variant calls by each amplification method at sites of low complexity (mainly TATA repeats). We also assessed the ability of each amplification condition to find a set of known consensus errors in the reference, and used data obtained to score true and false positives. The known consensus errors in the reference genome were previously identified and published [16]. We used the iCORN program [16] to conservatively map reads to the reference genome. iCORN iteratively maps reads and compares them with perfectly mapping reads. PCR-free variant calls were included as a quality control. A mismatch or indel was called with mpileup (a function in SAMtools) at a quality ≥40 from the vcf file. We differentiate between a true positive mismatch or indel (if the call agrees with iCORN and PCR-free), and a false positive, (if the call disagrees with iCORN and PCR-free) at that position. Using the total number of errors/variants counted per sample, we generated a normalized error score for each enzyme or library-preparation method.

Rank

We ranked the results by assigning a score of 1 to the best and a score of 0 to the worst. For each amplification condition, we calculated the average score value obtained in both coverage and fidelity metrics. Conditions with the highest average score was ranked highest while that with the least assigned the lowest rank (see Table 2).

Chimeric and duplicate reads analysis

We used the BWA paired-end alignment algorithm to generate mapped data sets. Picard software http://picard.sourceforge.net/command-line-overview.shtml was used to identify/mark duplicates from the mapped data set. Chimeric reads were identified from the alignment results using the “bwasw” option with single-end reads. Each replicate data set was analyzed independently to determine the number of chimeric or duplicate reads. We calculated the number of chimeric or duplicate reads as a percentage of the total mapped reads for each data set. Mean values are shown (see additional file 1, Table S2).

All datasets used in this study have been deposited in the ENA read archive under accession number [ERP000832].