Swab-Seq: A high-throughput platform for massively scaled up SARS-CoV-2 testing

Sample Collection

All patient samples used in our study were deidentified. All samples were obtained with UCLA IRB approval. Nasopharyngeal samples were collected by health care providers from individuals whom physicians suspected to have COVID19.

Creation of Contrived Specimens

For the clinical limit of detection experiments, we pooled confirmed, COVID-19 negative remnant nasopharyngeal swab specimens collected by the UCLA Clinical Microbiology Laboratory. Pooled clinical samples were then spiked with ATCC Inactivated Virus (ATCC 1986-HK) at specified concentrations and extracted as described below. For the clinical purified RNA samples, they were collected as nasopharyngeal swabs and purified using the KingFisherFlex (Thermofisher Scientific) instrument using the MagMax bead extraction. All extractions were performed according to manufacturer’s protocols. For extraction-free samples, we first contrived samples at specified concentrations into pooled, confirmed negative clinical samples and diluted samples in TE buffer or water prior to adding to the RT-PCR master mix.

Processing of Extraction-Free Saliva Specimens

Direct saliva is collected into a Matrix tube (Thermofisher, 3741-BR) using a small funnel (TWDRer 6565). The saliva samples were collected into a matrix tube and heated to 95°C for 30 minutes. Samples were then either frozen at −80C or processed by dilution with 2X TBE with 1% Tween-20, for final concentration of 1x TBE and 0.5% Tween-20²². We also tested 1x Tween with Qiagen Protease and RNA Secure (ThermoFisher), which also works but resulted in more sample-to-sample variability and required additional incubation steps.

Processing of Extraction-Free Nasal Swab Lysates

All extraction-free lysates were inactivated using a heat inactivation at 56°C for 30 minutes. Samples were then diluted with water at a ratio of 1:4 and directly added to mastermix. Dilution amounts varied depending on the liquid media that was used. We found that of the CDC recommended media, normal saline performed the most robustly. Viral Transport Media and Amies Buffer showed significant PCR inhibition that was difficult to overcome, even with dilution in water. We recommend placing the swab directly into the diluted TE buffer, which has little PCR inhibition.

Barcode Primer Design

Barcode primers were chosen from a set of 1,536 unique 10bp i5 barcodes and a set of 1,536 unique 10bp i7 barcodes. These 10 bp barcodes satisfied the criteria that there is a minimum Levenshtein²⁷ distance of 3 between any two indices (within the i5 and i7 sets) and that the barcodes contain no homopolymer repeats greater than 2 nucleotides. Additionally, barcodes were chosen to minimize homo- and hetero-dimerization using helper functions in the python API to Primer3²⁸. Additional details and code for primer design can be found at https://github.com/octantbio/SwabSeq.

Construction of S2 and RPP30 in vitro standard

RT-PCR was performed using primers shown above on gRNA of SARS-CoV-2 (Twist BioSciences, #1) for construction of a in vitro S2 standard DNA template. RT PCR (FP_1, R) and a second round of PCR (FP_2, R) was performed on HEK293T lysate for construction of a in vitro RPP30 standard DNA template. Products were run on a gel to identify specific products at ~150 bp. DNA was purified using Ampure beads (Axygen) using a 1.8 ratio of beads:sample volume. The mixture was vortexed and incubated for 5 minutes at room temperature. A magnet was used to bind beads for 1 minute, washed twice with 70% EtOH, beads were air-dried for 5 minutes, and then removed from the magnet and eluted in 100 μL of IDTE Buffer. The bead solution was placed back on the magnet and the eluate was removed after 1 minute. DNA was quantified by nanodrop (Denovix).

This prepared DNA template was used for standard HiScribe T7 in vitro transcription (NEB). IVT reactions prepared according to the manufacturer’s instructions using 300 ng of template DNA per 20 uL reaction with a 16 hour incubation at 37°C. IVT reactions were treated with DNAseI according to the manufacturer’s instructions. RNA was purified with an RNA Clean & Concentrator-25 kit (Zymo Research) according to the manufacturer’s instructions and eluted into water. RNA standard was quantified both by nanodrop and with a RNA screen tape kit for the TapeStation according to the manufacturer’s instructions (Agilent) to verify the RNA was the correct size (~133 nt).

One-Step RT-PCR

RT-PCR were performed using either the Luna® Universal One-Step RT-qPCR Kit (New England BioSciences E3005) or the TaqPath™ 1-Step RT-qPCR Master Mix (Thermofisher Scientific, A15300) with a reaction volume of 20 μL. Both kits were used according to the manufacturer’s protocol. The final concentration of primers in our mastermix was 50 nM for RPP30 F and R primers and 400 nM for S2 F and R primers. Synthetic S2 RNA was added directly to the mastermix at a copy number of 500 copies per reaction. Sample was loaded into a 20 μL reaction. All reactions were run on a 96- or 384-well format and thermocycler conditions were run according to the manufacturer’s protocol. We observe significant differences between the amplification of samples from purified RNA versus extraction-free (unpurified swab) samples (Figure S12). For purified RNA samples we performed 40 cycles of PCR. For extraction-free samples, we performed endpoint PCR for 50 cycles.

Multiplex Library Preparation

After the RT-PCR reaction, samples were pooled using a multichannel pipet or Integra Viaflow Benchtop liquid handler. 6 μL from each well were combined in a sterile reservoir and transferred into a 15 mL conical tube and vortexed. 100 uL of the pool was transferred to a 1.7 mL eppendorf tube for a double-sided SPRI cleanup ²⁹. Briefly, 50 μL of AmpureXP beads (Beckman Coulter A63880) were added to 100 μL of the pooled PCR volume and vortexed. After 5 minutes, a magnet was used to collect beads for 1 minute and supernatant transferred to a new eppendorf tube. An additional 130 μL of Ampure XP beads were added to the 150 μL of supernatant and vortexed. After an additional 5 minutes, the magnet was used to collect beads for 1 minute and the beads were washed twice with fresh 70% EtOH. DNA was eluted off the beads in 40 μL of Qiagen EB buffer. The magnet was used to collect beads for 1 minute and 33 μL of supernatant was transferred to a new tube. Purified RNA was quantified and library quality was assessed using the Agilent TapeStation. We observe some differences in non-specific peaks in our TapeStation analysis of the final library preparation, particularly when sequencing unpurified samples out to 50 cycles (Figure S13). The presence of non-specific reads affects the quantification of the library, loading concentration, and cluster density. Therefore, we suggest quantifying the final library based on the proportion of the desired peaks.

Sequencing Protocol

Libraries were sequenced on either an Illumina MiSeq (2012) or Nextseq 550. Prior to each MiSeq run, a bleach wash was performed using a sodium hypochlorite solution (Sigma Aldrich, 239305) according to Illumina protocols. We also perform a maintenance wash between each run. The pooled and quantitated library was diluted to a concentration of 6 nM (based on Qubit 4 Fluorometer and Illumina’s formula for conversion between ng/μl and nM) and was loaded on the sequencer at either 25 pM (MiSeq) or 1.5 pM (NextSeq). PhiX Control v3 (Illumina, FC-110-3001) was spiked into the library at an estimated 30-40% of the library. PhiX provides additional sequence diversity to Read 1, which assists with template registration and improves run and base quality.

For this application, the MiSeq requires 2 custom sequencing primer mixes, the Read1 primer mix and the i7 primer mix. Both mixes have a final concentration of 20 μM of primers (10 μM of each amplicon’s sequencing primer). The NextSeq requires an additional sequencing primer mix, the i5 primer mix, which also has a final concentration of 20 μM. The MiSeq Reagent Kit v3 (150-cycle; MS-102-3001) is loaded with 30 μL of Read1 sequencing primer mix into reservoir 12 and 30 μL of the i7 sequencing primer primer mix into reservoir 13. The NextSeq 500/550 Mid Output Kit is loaded with 52 μl of Read1 sequencing primer mix into reservoir 20, 85 μl of i7 sequencing primer mix into reservoir 22, and 85 μl of i5 sequencing primer mix into reservoir 22. Index 1 and 2 are each 10 bp, and Read 1 is 26 bp.

Analysis

The bioinformatic analysis consists of standard conversion of BCL files into FASTQ sequencing files using Illumina’s bcl2fastq software (v2.20.0.422). Demultiplexing and read counts per sample are performed using our custom software. Here read1 is matched to one of the three expected amplicons allowing for the possibility of a single nucleotide error in the amplicon sequence. The hamming distance is the number of positions at which the corresponding sequences are different from each other and is a commonly used measure of distance between sequences. Samples are demultiplexed using the two index reads in order to identify which sample the read originated from. Observed index reads are matched to the expected index sequences allowing for the possibility of a single nucleotide error in one or both of the index sequences. The set of three reads are discarded if both index1 and index2 have hamming distances greater than 1 from the expected index sequences. The count of reads for each amplicon and each sample is calculated. In this analysis we make use of a few custom scripts written in R that rely on the ShortRead ³⁰ and stringdist ³¹ packages for processing fastq files and calculating hamming distances between observed and expected amplicons and indices. This approach was conservative and gave us very low level control of the sequencing analysis. However, we anticipate that continued development of the kallisto and bustools SwabSeq analysis tools ¹⁹ will be a more user-friendly and computationally efficient solution for SwabSeq.

Criteria for Classification of Purified Patient Samples

For our analytic pipeline, we developed QC metrics for each type of specimen. For purified RNA, we require each sample to have at least 10 reads detected for RPP30 and that the sum of S2 and in vitro S2 standard reads exceeds 2,000 reads. If these conditions are not met, the sample is rerun one time and if there is a second fail we request a resample. To determine if SARS-CoV-2 is present, we calculate if the ratio of S2 to in vitro S2 standard exceeds 0.003. (We note that we add 1 count to both S2 and in vitro S2 standard before calculating this ratio to facilitate plotting the results on a logarithmic scale.) If the ratio is greater than 0.003 we concluded that SARS-CoV- 2 is detected for that sample and if it is less than or equal to 0.003 we conclude that SARS-CoV-2 is not detected (Figure S5C).

The same pair of primers will amplify both the S2 and in vitro S2 standard amplicons. Because we run an endpoint assay, the primers will be the limiting reagent to continued amplification. In developing this assay, we observed that as S2 counts increase for a sample, the in vitro S2 standard counts decrease (Figure S3). We found that at very high viral levels, in vitro S2 standard read counts decreased to less than 1000 reads. Therefore, analysis of S2 and in vitro S2 standard together allowed our QC to call SARS-CoV-2 even at extremely high viral levels.

Since the S2 and in vitro S2 standard are derived from the same primer pair, to account for the scenario where in vitro S2 standard counts are low because S2 amplicon counts are very high and the sample contains large amounts of SARS-CoV-2 RNA (Figure S3) in the QC we require that the sum of S2 and in vitro S2 standard counts together exceeds 2000. For example, if we detected greater than 2000 S2 counts and 0 in vitro S2 standard counts this would certainly be a SARS-CoV-2 positive sample and we would result: SARS-CoV-2 detected.

Criteria for Classification of Purified Patient Samples

We require that the sum of S2 and S2 synthetic spike-in reads exceeds 500 reads, or the results are considered inconclusive. With inhibitory lysates we have observed that this lower total is acceptable for maintaining sensitivity and specificity while limiting the number of tests that are considered inconclusive. If the sum of S2 and S2 synthetic spike-in reads exceeds 500, we determine if SARS-CoV-2 is detected in a sample by seeing if the ratio of S2 to S2 spike exceeds 0.003. (We note that we add 1 count to both S2 and S2 spike before calculating this ratio to facilitate plotting the results on a logarithmic scale.) If the ratio is greater than 0.003 we concluded that SARS-CoV-2 is detected for that sample and if it is less than or equal to 0.003 we conclude that SARS-CoV-2 is not detected as long as 10 reads are detected for RPP30 for that sample. This serves as a control that sample collection took place properly and contains a human specimen. If fewer than 10 reads are detected for RPP30 and the ratio of S2 to S2 spike is less than or equal to 0.003 the results are considered inconclusive. We have modified this criteria such that only samples without SARS-CoV-2 signal are considered inconclusive if RPP30 reads are fewer than 10. This ensures that we do not miss SARS-COV-2 positive samples that may have fewer RPP30 reads (Figure S7).

Downsampling analysis

Reads were downsampled from the results for the NP purified confirmatory LoD shown in Figure S5B. We observed that downsampling down to 5,000 reads per well resulted in no instances of mis-classification of SARS-CoV-2 presence or absence. At 5,000 reads per well approximately 3% of wells would no longer pass the filter that the sum of S2 and S2 spike reads exceeded 1,000 reads and would result in a sample being classified as ‘Inconclusive’. A logistic regression classifier described elsewhere¹⁹ should robustly tolerate a small fraction of outlier samples with slightly lower read depth.

Analysis of index mis-assignment

Unique dual indices and amplicon specific indices were used to study index mis-assignment. In this scheme, each sample was assigned two unique indices for the S2 or Spike amplicon and two unique indices for the RPP30 amplicon for a total of four unique indices per sample. A count matrix with all possible pairwise combinations, i.e. a “matching matrix”, was generated for each index pair (one i7 and one i5) using kallisto and bustools ³². The counts on the diagonal of the matching matrix correspond to input samples and counts off of the diagonal correspond to index swapping events. The extent of index mis-assignment for the i7 and i5 index was determined by computing the row and column sums, respectively, of the off-diagonal elements of the matching matrix. The observed rate of index swapping to wells with a known zero amount of viral RNA was determined by computing the mean of the viral S counts to spike ratio for those wells.

Improving Limit of Detection Requires Minimizing Sources of Noise

One of the major challenges in running a highly sensitive molecular diagnostic assay is that even a single contaminant or source of noise can decrease the test’s analytical sensitivity. In the process of developing SwabSeq, we observed S2 reads from control samples in which no SARS-CoV-2 RNA was present (Figure 1D). We subsequently refer to these reads as "no template control” (NTC) reads. A key part of SwabSeq optimization has been understanding and minimizing the sources of NTC reads in order to improve the limit of detection (LoD) of the assay. We identified two important sources of NTC reads: molecular contamination and mis-assignment of sequencing reads.

To minimize molecular contamination, we followed protocols and procedures that are commonly used in molecular genetic diagnostic laboratories³³. To limit molecular contamination, we use a dedicated hood for making dilutions of the synthetic RNA controls and master mix. At the start of each new run, we sterilize the pipettes, dilution solutions, and PCR plates with 10% bleach, followed by UV-light treatment for 15 minutes.

To prevent post-PCR products that are at high concentration from contaminating our pre-PCR processes, we physically separated pre- and post-PCR steps of our protocol into two separate rooms, where any post-PCR plates were never opened within the pre-PCR laboratory space. To further protect from post-PCR contamination, we compared RT-PCR mastermixes with or without Uracil-N-glycosylase (UNG). The presence of UNG in the TaqPath™ 1-Step RT-qPCR Master Mix (Thermofisher Scientific) showed a significant improvement reducing post-PCR contamination of S2 reads present in the negative patient samples as compared with the Luna One Step RT-PCR Mix (New England Biosciences) (Figure S14). The RT-PCR mastermix contains a mix of dTTP and dUTP such that post-PCR amplicons are uracil containing DNA. These post-PCR that are remnants of previously run SwabSeq experiments therefore can be selectively eliminated by UNG. Importantly, this addition does not interfere with downstream sequencing.

A third source of molecular contamination was carryover contamination on the sequencer template line of the Illumina MiSeq³⁴. Without a bleach maintenance wash, we found that indices from the previous sequencing run were identified in a subsequent experiment where those indices were not included. While the number of reads for some indices were present at a number of S2 reads, the presence of carryover contamination affects the sensitivity and specificity of our assay. After an extra maintenance and bleach wash, we substantially reduced the amount of carryover reads present to less than 10 reads (Figure S15).

Another source of NTC reads is mis-assignment of amplicons. Mis-assignment of amplicons occurs when sequencing (and perhaps at a lower rate, oligo synthesis) errors result in an amplicon sequence that originates from the in vitro S2 standard but is mistakenly assigned to the S2 sequence within a given sample. Only 6 bp distinguishes S2 from in vitro S2 standard at the beginning of read 1. Sequencing errors can result in in vitro S2 standard reads being misclassified as S2 reads as error rates appear to be higher in the beginning of the read (Figure S16A). If computational error correction of the amplicon reads is too tolerant, these reads may be inadvertently counted to the wrong category. To reduce this source of S2 read misassignment, we use a more conservative thresholding on edit distance (Figure S16B). Future redesigns or extensions to additional viral amplicons should consider engineering longer regions of sequence diversity here.

An additional source of NTC reads is when S2 amplicon reads are mis-assigned to the wrong sample based on the indexing strategy. In our assay, individual samples are identified by pairs of index reads (Figure 1B). Mis-assignment of samples to the wrong index could occur if there is contamination of index primer sequences, synthesis errors in the index sequence, sequencing errors in the index sequences, or “index hopping” ³⁵.

We leveraged multiple indexing strategies in our development of SwabSeq, from fully combinatorial indexing (where each possible combination of i5 and i7 indices was used to tag samples in the assay) to unique-dual indexing (UDI) where each sample has distinct and unrelated i7 and i5 indices (Figure S14). However, the ability to scale can be limited due to the substantial upfront cost of developing that many unique primers. Fully combinatorial indexing approaches significantly expand the number of unique primer combinations. We have also explored a compromise strategy between fully combinatorial indexing and UDI where sets of indices are only shared between small subsets of samples. Such designs reduce the effect of sample mis-assignment while facilitating scaling to tens of thousands of patient samples (Figure S17). With a fully combinatorial indexing (Figure S17) we observed that NTC read depth was correlated with the total number of S2 reads summed across all samples that shared the same i7 sequence (Figure S18A). This is consistent with the effect of index hopping from samples with high S2 viral reads to samples that share the same indices. It is possible to computationally correct for this effect, for example using a linear mixed model (Figures S18B).

Finally, the challenges associated with combinatorial and semi-combinatorial indexing strategies can be mitigated by using unique dual indexing (UDI), a known strategy to reduce the number of index-hopped reads by two orders of magnitude³⁶. We have observed consistently lower S2 viral reads for negative control samples with this strategy. It also enables quantification of index mis-assignment by counting reads for index combinations that should not occur in our assay (Figure S19 A and B). The number of index hopping events is correlated with the total number of S2 + S2 spike reads (Figure S19 C and D), indicating that hopped reads are more likely to come from wells where the expected index has strong viral signal. We quantify the overall rate of hopping as 1-2% on a MiSeq, and expect this rate may be higher on patterned flow cell instruments.

There are many sources of noise in amplicon-based sequencing, from environmental contamination in the RT-PCR and sequencing steps to misassignment of reads based on computational correction and “index-hopping” on the Illumina flow cells. Preventing and correcting these sources of error considerably improves the limit of detection of the SwabSeq assay.

Supplementary Documents

Optimized Protocol.docx

Equipment List.xlsx

Sample Experimental Design Setup.xlsx

Supplementary Figures

Supplementary Figures.pdf