Detecting SARS-CoV-2 and its variant strains with a full genome tiling array

Sample collection

Eight deidentified SARS-CoV-2–positive clinical samples from the Wyoming Public Health Laboratory were obtained and used in this study. All samples were collected before May 2020. They are referred to as S1 to S8 in this study. No consent was necessary because only deidentified viral materials were used for this study; no human samples were collected or used.

Sample preparation for Illumina sequencing

Samples were prepared as previously described using the ARTIC sequencing methods. In brief, cDNA was prepared from total RNA extracted from clinical samples using SuperScript IV (SSIV, Thermo Scientific) and random hexamer priming. The resultant cDNA was amplified in two PCR reactions using the ARTIC Pool1 and Pool2 SARS-CoV-2 v3 primer sets and Q5 High-Fidelity DNA Polymerase (NEB). Following PCR, samples were purified using AMPure XP SPRI beads (Beckman Coulter). Illumina adaptors were added using the NEBNext® Ultra™ II DNA Library Prep Kit (NEB) and SPRI bead purification was repeated.

Illumina sequencing data processing

Sequencing data of SARS-CoV-2 were aligned using BWA 0.7.17 [11] using reference genome NC_045512 downloaded from NCBI. Binary alignment map (BAM) files were sorted and indexed using Samtools 1.9 [12]. BCFTools 1.9 was used to count allele frequency from the BAM files.

Sample preparation for chip hybridization

To prepare samples for hybridization to the chips, 0.05 μL of purified PCR product was reamplified using the ARTIC protocol and Pool1 and Pool2 v3 primer sets for 35 cycles with 50 μM biotin-11-dUTP (Jena Biosciences) added to the reaction mixture. Pool1 and Pool2 were combined for each sample and fragmented using DNase I (D4263, Sigma). About 2000 Kunitz units of lyophilized DNase I was resuspended on ice using 2 mL of 1x DNase I Buffer (10 mM Tris–HCl pH 7.5, 2.5 mM MgCl₂, 0.1 mM CaCl₂). The resuspended enzyme was diluted 1000-fold using 1x DNase I Buffer and an equal volume was added to samples prewarmed to 37°C. Samples were incubated for 30 min at 37°C and the reactions were stopped by adding EDTA to a final concentration of 12.5 mM and incubating for 20 min at 75°C.

Hybridization

About 45 μL of the fragmented sample was hybridized overnight at 45°C to the chip in a 60 μL final volume containing 5 mM EDTA, 6.25 mM HEPES pH 8.0, 312.5 mM NaCl, 1.25% FIcoll 400, 0.5 nM Cy3-AM1 (GCTGTATCGGCTGAATCGTA). Following hybridization, chips were washed for 10 min at room temperature in Wash A (2x SSC, 0.1% TWEEN-20) and then for 10 min at 39°C in Wash B (0.5xSSC, 0.1% TWEEN-20). Chips were stained for 15 min at room temperature using 0.02 mg/mL Cy3-Streptavidin (Thermo) in 4x SSC and washed for 5 min at room temperature using 4xSSC. Chips were scanned using a custom-built confocal scanner for 0.5, 1, 4 and 8 s in the green (Cy3) channel in 4x SSC.

SARS-CoV-2 tiling array design

A tiling array is based on traditional microarray technology but aimed to resequence the entire genome. Our tiling array was designed based on the SARS-CoV-2 genome (NC_045512), which contains ~30 K nucleotides. Each position in the SARS-CoV-2 genome is captured by eight probes, four for sense and four for antisense strands. The four probes for the sense or antisense can be considered as a probe set, compactly designated as a ‘probeset’ hereafter. A probe is a 25-mer synthetic oligonucleotide matching the SARS-CoV-2 genome. In each probeset, the only difference among the probes is the middle nucleotide, which is used to capture the four possible variants (A, T, C, G) at that specific position (Figure 1A). Sense and antisense probesets are designed based on the same concept except that the antisense probes are used to capture the same information from the antisense strand resulting from reverse transcribed double-strand cDNA. Because the probe design requires 12 flanking sequences from either side of the targeted base, the first and last 12 bases of SARS-CoV-2 genome are not covered by the tiling array. In total, there are 239 000 probes on the tiling array. Each probe occupies a 3μm² area and the entire tiling array is 3mm². Due to variations in the intensity within the array, we elected to scan the array with three different exposure times (0.5, 1 and 4 s). Longer exposure time allows more accurate intensity profiling on weak intensity regions. However, it also saturates the brighter regions on the array. A proper summation method is required to combine the results from the six replicates (two strands x three exposures) for a same genomic position.

Base calling algorithms

Considering the strandness and the exposure parameter, each nucleotide position in the SARS-CoV-2 genome is represented with six probesets (two strands × three exposures) on the tiling array. Within a probeset, the sorted intensity values returned by individual probes are represented as . The highest intensity identifies a probe of a specific nucleotide base, one from , and that specific nucleotide base is denoted as a tentative call from the probeset in question. To make an ultimate base call for each position in the genome, we developed two base calling methods: weighted voting and ranked credibility score (Figure 1B).

Evaluation metrics with running genome windows

We calculated three metrics for a segment of nucleotide sequence along the SARS-CoV-2 genome. The GC content is dependent on the reference genome of SARS-CoV-2. The sequencing depth is dependent on the high-throughput sequencing result. The replicate consistency is defined specifically for the SARS-CoV-2 tiling array data. GC content is the percentage of guanine (G) and cytosine (C) bases within the nucleotide segment. Sequencing depth is the average reads coverage within the nucleotide segment. Across all segments of the whole genome, the sequencing depth values were standardized to a minimum of 0 and a maximum of 1. Replicate consistency is first defined for one exact nucleotide position as n/6 with n being the number of replicate probesets supporting the ultimate position base call; furthermore, replicate consistency of a nucleotide segment is obtained as the average number across all the individual positions within the segment. Given the linear single strand of SARS-CoV-2 genome, we typically binned the genome in running, nonoverlapping windows of 100 bp wide, unless specified otherwise. The aforementioned metrics, GC content, sequencing depth and replicate consistency, are calculated for each running bin and are thus connected into trend lines along the whole genome.

Unanimity and consistency of base calls

The SARS-CoV-2 tiling array builds in six replicated intensity readings for each SARS-CoV-2 sample. For the eight samples we tested, the intraclass correlation coefficient [13] between the six replicate series ranged over 0.88–0.90 with an average of 0.89, and 79.8–87.1% (mean: 84.6%) of all 30 K genomic positions showed perfect replicate consistency, i.e. unanimous tentative calls among six replicates. For all genomic positions represented on a tiling array, we divided them into a probeset-unanimous group and a nonunanimous group. We calculated a noise statistics for each genomic position, which was defined as the reciprocal of relative highest intensity (see METHODS). Because SARS-CoV-2 is a haploid genome, theoretically, only one allele should be observed at a genomic position. Hence, within a probeset, we may assume only one probe represents the correct base, while the other three, summarized in the noise statistics, represent background noise. Aggregating all eight samples, we plotted the distribution of the noise statistics for the unanimous positions and the nonunanimous positions separately. To our expectation, genomic positions of unanimous tentative calls from replicate probesets demonstrate a strikingly lower noise level than positions of nonunanimous tentative calls (Figure 2A). In other words, unanimous genomic positions generally show low noise level, or striking predominant intensity level, within each probeset.

In the landscape of replicate consistency along the whole genome, several notable genomic regions exhibited minor consistency drops (Figure 2B). One such region was located between position number 19 300 and number 19 550, where distinctively lower consistency was visually recognizable from the preceding and succeeding levels. Interestingly, the sequencing depth trend line displayed evident dips at roughly the same genomic regions (Figure 2C). After averaging over the eight samples, the replicate consistency and sequencing depth had a Pearson correlation coefficient of 0.32 (P < 0.0001). This result suggests that tiling array hybridization and high-throughput sequencing were faced with the same difficulty along the SARS-CoV-2 genome. It is well known that GC content is associated with sequencing depth [14]. We summarized GC content of SARS-CoV-2 genome in sliding windows of 100-bp width. The overall GC content of SARS-CoV-2 genome is 38.0%, oscillating between 20% and 56% in sliding windows. Using linear regression, we found that GC content can explain a small portion of the variation in sequencing depth (adjusted R² = 0.036, P = 0.0006). Similarly, GC content was also associated with array consistency (adjusted R² = 0.13, P < 0.0001). It seems that GC content partially explains the genome region difficulties encountered by both tiling array hybridization and high-throughput sequencing. However, there are yet other undiscovered factors affecting the resolution of particular regions of the SARS-CoV-2 genome.

Validation via high-throughput sequencing

To validate the performance of SARS-CoV-2 tiling array, we conducted Illumina high-throughput sequencing on all eight samples. Using sequencing results as the gold standard, the tiling array displayed remarkable accuracy. On average, the eight samples had an accuracy of 99.61% (range: 99.50–99.78%) for the weighted voting method and 99.57% (range: 99.35–99.81%) for the ranked credibility score method, respectively (Figure 3). Such high accuracy certifies the high reliability of SARS-CoV-2 tiling array coupled with the base calling methods. The two base calling methods performed very similarly. The weighted voting method had a slight advantage in six of the eight samples, while the ranked credibility score method outperformed in the rest of the two samples. The two methods make the final decision via different mechanisms: one by voting and the other by ranking, and they each may make correct calls in different scenarios. For example, in sample S2, at position number 22 908, the weighted voting method voted for nucleobase A based on the weighted decision-making vector of (1.073-A, 1.035-A, 1.027-A, 1.049-A, 1.0532-T and 1.060-A); the ranked credibility score method called nucleobase T based on the credibility score vector of (0.507-A, 0.156-A, 0.669-A, 0.307-A, 0.938-T and 0.916-A). Sequencing result approved the weighted voting method for this genomic position. In another example, in sample S1, at position number 25 013, the weighted voting method voted for nucleobase C based on the weighted decision-making vector of (4.698-C, 2.768-G, 6.539-C, 3.123-G, 9.275-C and 3.782-G); the ranked credibility score method called nucleobase G based on the credibility score vector of (0.99993-C, 0.9994-G, 0.9988-C, 0.9988-G, 0.9995-C and 0.99995-G). Sequencing approved the ranked credibility score method for this genomic position.

Mutation detection by tiling array

Variant SARS-CoV-2 strains add a paramount danger to the ongoing SARS-CoV-2 pandemic. When the base call for a genomic position differs from the nucleobase specified in the reference genome, a point mutation is suggested. Of the eight samples, SARS-CoV-2 tiling array detected 61 mutations at 35 genomic positions. All of these 61 mutations were recovered in the sequencing results, indicating a sensitivity of 93.8% and a specificity of 100% (Figure 4A). On average, each sample showed eight mutations. SARS-CoV-2 has four proteins: spike, envelope, membrane and nucleocapsid. From an epidemiology perspective, the spike protein is the most important protein because it is the target of current vaccines. Our analysis identified two mutated positions in the coding region of the spike protein: position 23 403 was mutated in four samples, coincident with the infamous D614G mutation; position 24 694 was mutated in two samples, without changing the amino acid outcome (Figure 4B). Another protein, nucleocapsid, has increasing importance because it has been suggested as an alternative vaccine target [15]. Our analysis identified four mutations in the coding region for nucleocapsid protein, two nonsynonymous and two synonymous (Figure 4C). The first nonsynonymous mutation appeared in sample S6 at position number 28 311, changing Pro to Leu (P13L). The second nonsynonymous mutation also appeared in sample S6 but was located at position number 28 863, changing Ser to Leu (S197L).