In this study we tested the in silico performance of 43 PCR assays (Supplementary Table 1) with different gene targets [19 ORF1ab, 8 spike (S) protein, 1 ORF3a, 3 envelope (E) protein, 1 ORF8, and 11 nucleocapsid (N)] against a set of 1,690,689 SARS-CoV-2 accessions (Supplementary Table 2) and sequences downloaded from the GISAID EpiCoV™ database on July 7, 2021 (11, 13, 16–26). At the time of analysis, the CDC listed Pango lineages P.1, B.1.351, B.1.1.7, B.1.427, and B.1.429 as VoCs and lineages P.2, B.1.525, and B.1.526 as VoIs (27, 28). Additional follow-up analysis was performed on the B.1.1.529 VoC as it emerged later. The PSET definition for each assay target was based on a reference amplicon sequence with delimited primer and probe regions (Figure 1). Twenty nucleotides of additional sequence context outside the amplicon were also included to inform alignment at the 5′ and 3′ ends. Context and inter-primer regions were obtained via global-local alignment of the assay primers to the SARS-CoV-2 reference genome by running the glsearch36 program of the FASTA suite (29).

The first phase of PSET analysis queries the assay target definition against a BLAST+ database to search for matching subjects based on local alignment with the blastn program (30). Expansion of ambiguous DNA codes of the query is required for compatibility. For example, “GAWTAYA” has two, representing four expansions: “GAATACA,” “GATTACA,” “GAATATA,” and “GATTATA.” Only the first permutation is queried. An additional step re-evaluates the BLAST+ identity statistics by replacing the expanded query with the original one to account for ambiguous base similarity. Subjects with ≥85% identity to the amplicon region are then extended to cover the query range and extracted.

The second phase queries the corresponding primers separately against the library of extracted sequences to search for matching subjects based on global-local alignment with the glsearch36 program, which is compatible with ambiguous DNA codes. Subjects with ≥90% identity are kept and aggregated by subject accession. A true positive (TP) is called if all primers aligned to the subject with the required identity and strand arrangement such that the primer would hypothetically amplify the target. The special case of a perfect true positive (PT) is called when there is 100% identity. Otherwise, a false negative (FN) is called. If the subject contains an N in any of the primer alignment regions, an N is appended to the category name (TPN and FNN). Doing so avoids removing low-quality, yet potentially informative sequences. An additional unknown (UNK) category is included to account for subjects that failed alignment during the first phase. In the case of near-neighbor analysis, a false positive (FP) or true negative (TN) is called when the taxonomy of the subject differs from the assay target. FPs with 100% identity are perfect false positives (PF).

Results were filtered to include high-coverage, human-host sequences with assigned Centers for Disease Control (CDC) and Pango lineages corresponding to VoC [Alpha (B.1.1.7), Beta (B.1.351), Gamma (P.1), Delta (B.1.617), Epsilon (B.1.427/B.1.429)] and VoI [Eta (B1.525), Theta (P.2), and Iota (B.1.526)] (27, 28). Metadata also included VoC and VoI designation times from disease control centers (31–39). Further aggregation and visualization was performed using the R programming language with the tidyverse (v1.3.1), lubridate (1.8.0), tsibble (v1.1.0), and cowplot (v1.1.1) packages (40–43). The PSET workflow itself was implemented using Biopython and Snakemake (44, 45). A summary of results and methodology refinements were posted to virological.org on a near-weekly basis (13).

For variant analysis, the local GISAID EpiCoV™ database used for the PSET analysis was filtered to include only unambiguous DNA codes, resulting in a subset of 961,051 sequences. Single nucleotide variations (SNVs) and insertions or deletions (indels) were calculated by running the nucmer and show-snps programs of the MUMmer4 suite (46). Parameters for nucmer included NCBI GenBank accession no. NC_045512.2 (47) as the reference and flags to search the forward strand of each query (-f) with 28 threads (-t). The show-snps program transformed the resulting delta file into a tabular file (-T) of SNP/indel calls. An R script loaded the variants identified via MUMmer4 and split the data into two groups according to overlap with assay primer and probe target regions. The percentage of positions with ≥n mutations was calculated, where n was in the range [0, 500]. A two-sample Kolmogorov-Smirnov procedure tested the null hypothesis that the observed mutation percentages arose from the same distribution for the assay target and non-target genomic regions.

The Omicron (B.1.1.529) wave emerged after the initial analysis was completed. Accordingly, an updated analysis targeted this VoC. On February 16, 2022, genomes were downloaded manually from the EpiCoV™ search page. Filters were enabled to include complete sequences with high coverage and collection dates. The PSET workflow analyzed the final set of 48,358 sequences (Supplementary Table 3) as previously described.

Figure 2 shows the target region of each assay with respect to the reference SARS-CoV-2 genome (NCBI GenBank accession no. NC_045512.2), SNVs, and indels (47). Mutations with <1% prevalence (~99.62% of all positions with observed mutations) were excluded from the figure but remain listed in Supplementary Table 4. The figure also plots the assay target regions on the genome. Plotting all the regions together with respect to mutations helps visualize signature erosion due to genetic drift.

On the other hand, when the assay target regions are compared specifically against the non-target regions, assay target regions with mutations are present in higher percentage of genomes (Figure 3). The two-sample Kolmogorov-Smirnov test obtained a p value of 2.2E−16, rejecting the null hypothesis at a significance level of 0.05 in favor of the two-sided alternative.

The assays were able to correctly detect the vast majority of the 1,690,689 genome sequences according to the in silico PSET analysis. Of the 43 assays, 34 aligned with 100% identity to over 1.6 million subjects. Otherwise, TP rates exceeded 98%, except for the Young-S and China_N assays, which exhibited high FN rates at ~47% and ~59% respectively. The TP percentage was calculated as the sum of PT, TP, and TPN divided by the total number of subjects. Alignments exceeded identity threshold for most corresponding subject amplicon sequences, only failing on average ~0.74% (UNK%) of the time. Table 1 and Supplementary Table 7 include a full breakdown by assay. In the separate Omicron study, most of the assays (40 / 43) were able to detect sequences with a TP rate of ≥96%, with Young_S, Thailand_WH-NIC_N, and China_N assays having reached TP rates of 75.64%, 2.98%, and 0.34% respectively (Supplementary Table 8). We also looked at the specificity of the assays against near-neighbor sequences (Supplementary Table 9), where 8 showed significant false positives (FP) and perfect FPs (PFs) (those with 100% identity).

Alignment identity of all assays against sequences of different lineages of SARS-CoV-2 was assessed from 2020-03-15 to 2021-07-05. Figure 4 depicts heat maps of the TP rate over time for each VoC/VoI with a corresponding line graph of the cumulative log-total number of sequences. Only GISAID sequences with collection date metadata specifying year, month, and day were included in the heatmap. The graph reveals sudden changes in assay alignment identity with respect to variant abundance.

Alignment identity remained constant for most assays, with perfect or near-perfect TP rates. However, some interesting patterns were observed. Some assays targeting the S and N genes performed poorly. For example, Young-S exhibited low TP rates within the B.1.1.7 and B.1.525 lineages and faltered for B.1.258, B.1.427, and B.1.429 isolates as cumulative sampling increased. Other assays targeting the S protein gene performed well, despite some temporary rough patches for Chan-S, Won-S and C5_COV_S_gene. Recently, Chan-S appears to have started failing for the B.1.258 lineage. China_N failed for the B.1.1.7, P.2, and P.1 clades. All other assays targeting the N gene performed well. However, Japan_NIID_2019-nCOV_N, HKU-N, cdc_n2, and Chan-N recently started failing for the B.1.525 lineage. Also, Young-N temporarily exhibited a low TP rate before recovering.

Figure 4 displays the timeline for each lineage to be identified as a VoC/VoI by the CDC, European CDC (ECDC), Public Health England (PHE), and World Health Organization (WHO). Most identifications occurred within a month of each other. However, for P.2 and B.1.429, 186 and 118 days passed between the initial and final call. The CDC and PHE called the former within 12 days of each other and the former within 111 days. In both cases, the WHO called last. We found no evidence of a call from the ECDC in these cases. The B.1.429 case shows the initial call occurring during a period of low alignment identity for the Young-S assay and the beginning of a decline for the Japan_NIID_WF-1_Seq_F519 assay. Despite the varying times at which these lineages were identified, when looking at the cumulative increase of the GISAID sequences for each, it is readily apparent that they initially occurred at very low proportions. This increase, and the resulting VoC/VoI designations, underscore the importance of continual monitoring from multiple organizations offering different perspectives and criteria with respect to spatiotemporal trends.

As assay targets were distributed throughout the genome, we were also interested in seeing whether emerging variants possessed mutations that resulted in more than one assay producing a predicted FN. Supplementary Table SB shows that nearly all (>99%) sequences resulted in a predicted true positive in at least 41 of 43 assays. Also, Supplementary Table SC shows that the majority (97%) of sequences producing a single FN were observed to originate from the China_N (87%) or Young-S (10%) assays. The majority (765,813) of sequences produced FNs in two different assays, China_N and Young-S, while 643,876 sequences did not produce an FN in any of the 43 assays. We observed a substantial drop-off in the number of sequences which caused FNs in 3 or more assays. This suggests that it is unlikely for a sequence to produce FNs in more than 2 assays and that a multiplex panel with any of the other 41 assays evaluated here would likely perform successfully.

The rapid increase of novel SARS-CoV-2 variants raises concerns regarding the efficacy of PCR-based diagnostic assays. This study assesses the potential signature erosion of the current COVID-19 real-time diagnostic assays using the GISAID sequence database. We find that 41 out of the 43 assays continue to perform very well throughout the observed timeframe even as new variants arose. Additional analysis completed after the emergence of Omicron revealed that 40 of 43 assays continued to perform well.

Variants containing mutations within the target regions for primer and probe hybridization were observed to have a substantial effect on predicted assay alignment identity. The characteristic S:Δ69/70 mutation of the Alpha variant occurs within the forward primer of the Young-S assay target, while mutations N:R203K and N:G204R occurred within the borders of the China_N assay target, resulting in PSET calling significant FN rates of those assays. The China_N assay failed completely for the lineages B.1.1.7, P1, and P2. This may be due to a “GGG” to “AAC” mutation in the forward primer, as observed in previous studies (8, 9). The 3-nt substitution mutation could potentially reduce the China_N forward primer binding affinity, which could significantly increase the FN rate, especially since it occurs at the 5′-end. The B.1.1.529 lineage contains an N:Δ31/33 deletion (48) which falls inside the borders of the Thailand_WH-NIC_N assay and results in PSET predicting significant FN rates. This may be a similar situation as described above where primer binding affinity is affected by this type of mutation.

Various point mutations at the primer binding motifs are observed in all the lineages, which did not affect overall performance according to evaluation of in silico alignments. The exception was for lineage B.1.258, which bears a complementary sequence to the China_N primer resulting in perfect alignment. The China_N and Young-S assays experience 56 × and 44 × the FN rate of all other assays combined, respectively. In almost all assays, there is a percentage of mutations in one of the targets of the assay.

The COVID-19 pandemic has demonstrated the need for the in silico monitoring of PCR-based diagnostic assay performance. Furthermore, given that surveillance network data are widely available to aid in quick-turnaround research and development of replacement assays, in silico monitoring has become an important early warning indicator (10). Our approach is consistent with recent alignment-based studies by confirming the overall high assay target sequence identity and detection of the China_N mutation (8, 9). We also observed potential early signs of signature erosion prior to official variant calls. Future research can further elucidate validation of in vitro and in silico predictions while considering the clinical and public health relevance. For example, the PSET algorithm could incorporate mismatch position, since mismatches near the 3′ end of the probe sequence potentially have a greater impact on assay performance (49, 50). Additional wet lab experiments can help systematically select alignment thresholds and motivate algorithmic refinements.