The variability and reproducibility of whole genome sequencing technology for detecting resistance to anti-tuberculous drugs

Samples

DNA was extracted and purified from the liquid cultures using a cetyltrimethylammonium bromide (CTAB) method [14]. The quality was assessed by fluorometric quantification, Qubit™ 3.0 Fluorometer with a dsDNA Broad Range Assay Kit (Thermo Fisher Scientific) and agarose gel electrophoresis. Triplicate DNA samples from each clinical isolate were prepared (biological replicates) and individual DNA extracts were subjected to repeated sequencing (technical replicates).

Library preparation and sequencing

For MiSeq sequencing, ~200 ng of genomic DNA was sheared to an average size of 500 bp by ultrasonication (Covaris S220). Sheared DNA was purified/concentrated on MinElute Spin Columns (Qiagen). DNA concentrations were measured on a Nanodrop UV spectrophotometer and the sheared samples diluted to 5–12.5 ng/μl. Library constructions were performed using the Ovation Rapid DR Multiplex System (NuGen) according to the manufacturer’s instructions. Purified libraries were amplified in emulsion PCR, size selected (500–700 bp) by preparative electrophoresis on composite gels (1.2% LMP-Agarose/0.8% Synergel) and then purified on MinElute Columns. Libraries were sequenced with an Illumina MiSeq V3 and 300-bp paired-end reads with samples randomised across two runs (each ~24 h in duration).

Ion Torrent library preparation and sequencing was performed at Thermo Fisher Scientific, UK. Sequencing was carried out with the Ion Torrent PGM™ system (Ion PGM). Libraries were constructed with the Ion Xpress™ Plus Fragment Library Kit as per the manufacturer’s instructions (MAN0009847 Revision C.0), using 100 ng of genomic DNA which was sheared with the provided Ion Shear™ Plus Reagents to an average size of 350 bp, size selection using an E-Gel® SizeSelect™ 2% Agarose Gel, and purification with Agencourt® AMPure® XP Reagent. Finally, the libraries were quantified on the StepOnePlus™ System using the Ion Library Quantitation Kit, then diluted to 100 pM and pooled in equal volume. Purified libraries were sequenced with an Ion 318™ v2 chip (400-bp kit) and the Ion PGM™ HiQ™ Chef Kit as stated in the manual (MAN0010919, revision A.0). The runtime was ~3 h per sample. The software used on both Ion PGM™ and the Ion Chef™ System was Torrent Suite™ Software version 4.6.

Bioinformatic pipeline

For the bioinformatic analysis we used a previously reported pipeline [10, 15, 16]. Unless stated otherwise, software was run at default settings. Reads were trimmed by Trimmomatic using a PHRED quality of 20 as the cutoff. Trimmed reads were then mapped against H37Rv (GCA_000195955.2) with BWA-mem v0.7.12 [17]. SNP and insertion and deletion (indel) variants were called with Samtools 0.1.19 [18] and GATK v3.6 [19]. We compared the variants called by both algorithms, but also report results of the conservative and typical approach of retaining the consensus polymorphisms across both methods. The genotypes of SNPs were called when an alternative allele was found in 20% of the mapped reads at a particular position. A default minimum depth of ten reads was required to call SNP genotypes, otherwise genotypes were denoted as missing data. This cutoff has been applied widely [15, 16, 20]. The robustness of the genotype calls was assessed across a range of depths of coverage of the reference and alternative alleles (depth 5–20, major allelic frequency >0.5 and >0.7). The reference genome was partitioned into overlapping 300-bp sequences allowing the uniqueness of genomic regions to be determined using gem-mappability [21]. Only 1.5% of the genome was estimated to be non-unique, and variants within these regions were discarded, leaving a set of high quality SNPs and indels. All 36 candidate drug-resistance genes [5] were found to be unique, thus removing the risk of false calling of SNPs due to inappropriate mapping to an analogous region. A summary of the pipeline is presented in Additional file 1: Figure S1.

In silico profiling of M. tuberculosis resistance phenotypes

We compared two informatics tools for assigning resistance from sequence data. Drug-resistance status across 14 drugs was called in silico from raw sequence data using the web-based TBProfiler tool (http://tbdr.lshtm.ac.uk/) [5]). This tool also generates lists of mutations in candidate loci, and these formed the basis of identifying any additional putative novel polymorphisms. All mutations were checked by analysis of alignments and de novo assembly, as well as confirmed by alternative sequencing methods (see the next section, “Confirmation of mutations detected by whole genome sequencing”). Resistance profiles were also generated with the Mykrobe predictor tool (version July 2016) [22].

Confirmation of mutations detected by whole genome sequencing

Genomic DNA was extracted as described above and used for PCR amplification prior to examination by line probe assay and/or DNA sequencing. The Genotype MTBDRplus (Hain Lifescience) investigates the rpoB and katG genes and inhA regulatory region and Genotype MTBDRsl (version 1, Hain Lifescience) investigates rrs, gyrA and embB. Both kits were used according to the manufacturer’s instructions. As the line probe assays encompass a limited number of loci, we also performed Sanger sequencing for inhA, katG, tlyA, eis, gidB, pncA, gyrA, ethA, embB, embC-embA, rpsL, folC and thyX genes (see Additional file 2: Table S1 for the primers used). PCR products were purified and both strands sequenced at StabVida (Portugal). All sequences were edited and analysed with ChromasPro 2.0.0 (Technelysium, Australia), compared to the sequences of M. tuberculosis H37Rv reference strain (GenBank AL123456.2) and aligned with Clustal Omega [23].

Coverage

Triplicate “extraction” DNA samples from ten clinical isolates and a single H37Rv sample were sequenced on the MiSeq platform. Four DNA samples (from POR5, 6 and 7 and H37Rv) were each sequenced six times (“technical” replicates). Duplicate DNA samples from three clinical isolates (POR1, 2 and 6) were also sequenced on the Ion PGM. Summaries of the sequence data obtained for each platform are presented in Additional file 3: Table S2. With MiSeq sequencing the number of paired reads varied across samples (median 1.2 million, range 0.4 to 3.2 million), and on average 99% of reads mapped to the H37Rv reference, giving a median depth of coverage of 51-fold (across sample range 18- to 79-fold). The majority of the genome (>96%) was covered to at least tenfold depth.

In contrast, for the Ion PGM platform the median number of reads was 990,854 (range 928,006–1,124,215) translating into a median of 53-fold (range 48- to 59-fold) genomic coverage. A large proportion of the genome (~25%) had low coverage and was attributed to regions with high GC content (Fig. 1). Whilst high coverage (100- to 200-fold) was attained for regions with GC content up to 69%, above this level coverage drops below tenfold, which was the cutoff used for calling variants. For MiSeq sequence data, this drop only occurs when the GC content reaches 75% or above. Many regions in the M. tuberculosis genome, especially the pe/ppe genes [24], are high in GC content (median 69%, range 47–87%) and therefore potentially difficult to characterise. The coverage across the 36 drug-resistance candidate genes was high for MiSeq (mean ~90-fold) and exceeded the tenfold cutoff, except in the thyA gene in the three POR1 replicates (Fig. 2). These XDR-TB replicates contained double dfrA-thyA deletions, thought to be responsible for para-aminosalicylic acid (PAS) resistance [25]. A direct comparison of the POR1, 2 and 6 sample coverage across platforms highlighted greater variability in candidate genes in Ion PGM due to differential GC content. Whilst there was platform-wide detection of the deletion-driven lower coverage in thyA in POR1 (Fig. 3; Additional file 4: Figure S2), the variable coverage in the neighbouring regions for Ion PGM could lead to less certainty in detection.

SNP variants and error rates

We estimated the variant error rates (measured as the number of sites which were discordant among replicates) to be low for both platforms (MiSeq 1 per 77 kbp, Ion PGM 1 per 41 kbp). Across comparable samples, the number of high quality SNPs detected using MiSeq data was higher than from Ion PGM, mostly due to low coverage in the alignments generated from the Ion PGM (Additional file 3: Table S2). We sought to investigate the effects of variant calling algorithms on the numbers of SNPs detected in unique genomic regions. From the MiSeq H37Rv data, similar numbers of SNPs were detected across replicates (Samtools 64–69 SNPs and GATK 69–79 SNPs, overlap 69 SNPs), supporting the existence of those variants and high sequence reproducibility (Additional file 5: Table S3). Across clinical isolate replicates the number of SNPs identified was similar and the overlap between variant calling algorithms was high (>90%; Additional file 5: Table S3). This observation was supported by the Ion PGM data but, due to uneven coverage, at least 120 SNPs fewer were identified when compared to matching MiSeq samples. Within platforms and calling algorithms there was variation between replicates in the indels detected, but there was high overlap between algorithms (>90%; Additional file 5: Table S3). Compared to SNPs the breakpoints for these variants are more difficult to characterise from alignments.

For the MiSeq platform data we assessed the number of SNP genotypes called across a range of coverage depths of the reference and alternative alleles (total depth 5- to 20-fold; major allelic frequency >0.5 and >0.7). The number of SNPs decreased pseudo-linearly with decreasing minimum read depth for H37Rv (87 to 67 SNPs; Additional file 6: Figure S3) and the ten clinical isolates (2290 to 2097 SNPs; Additional file 7: Figure S4). In general, differences in the number of SNPs between the Samtools and GATK algorithms decreased as the depth of coverage and allelic frequency thresholds increased. For H37Rv, read depths in excess of 20-fold had no impact on variants detected. Across clinical isolates, the highest possible stringency tested consisted of using a minimum coverage of 20 and an allelic frequency of 0.7 and led to near identical numbers of total SNPs called by the two variant calling algorithms (Samtools 1943, GATK 1990, either 2097, both 1840 SNPs; Additional file 7: Figure S4). Much of the discordance in the number of SNPs within replicate groups is due to differences in coverage leading to some polymorphisms not passing quality control filters. Using SNPs for which all replicates have non-missing genotypes, all replicates had identical numbers of SNPs except POR3C, which differed by two SNPs between POR3A and POR3B. Overall, the analyses indicated no major differences in SNPs detected between the two calling algorithms, and this supported the use of consensus variants for downstream analysis. For example, the set of common SNP variants were used to cluster all samples within a phylogenetic tree using FastTree v2.1.7 [26] (Additional file 8: Figure S5). Perfect clustering was observed amongst isolates and their replicates. At a finer resolution, we analysed the SNP differences between the replicates, and none were identified.

Calling in silico resistance phenotypes

When the MiSeq raw sequence data were subjected to analysis using TBProfiler, agreement with phenotypic susceptibility testing was high (95.3%, 82/86; Table 1). Discrepant results were recorded for PZA (×2) and PAS (×2) where phenotypically resistant isolates not identified by TBProfiler were found to have novel mutations in known candidate genes (Additional file 9: Table S4). The novel polymorphisms included a deletion in pncA of 20 bp (nucleotides 437–449) and a nucleotide insertion (GG) between codons 130 and 131. PAS-resistant isolates had a folC S98G mutation and a thyX G-4A, thyX I161T, dfrA-thyA deletion. Phenotypic testing of kanamycin drug susceptibility was not performed, but mutations associated with its resistance were detected in all eight isolates (Table 1; Additional file 9: Table S4). All mutations were confirmed using independent Sanger capillary sequencing and/or the line probe assays Genotype MTBDRplus and Genotype MTBDRsl (Hain). Phenotypic resistance profiles were confirmed and quantified by the qDST method for the MGIT960 system [12, 13].

The Mykrobe predictor tool was also applied to in silico call resistance. This approach looks for mutations associated with resistance to first-line drugs (rifampicin, isoniazid, ethambutol) and second-line drugs (streptomycin, ciprofloxacin, ofloxacin, moxifloxacin, amikacin, kanamycin, capreomycin). Of the 72 resistance calls made, 19 (26.4%) were incorrectly called “susceptible”. False negative calls were made for isoniazid (×1), ethambutol (×2), streptomycin (×4), amikacin (×4), and capreomycin (×3). Additionally there was a disagreement between TBprofiler and Mykrobe predictor with four samples for kanamycin, the latter program calling them as “susceptible” (Table 1).

For Ion PGM, when predicting individual drug-resistance profiles in the three isolates, in one isolate the gyrA D94A mutation associated with fluoroquinolone resistance could not be detected due to lack of coverage (Additional file 5: Table S3). However, the mutation was recovered if the coverage threshold was relaxed from ten- to fourfold.