SARS-CoV-2 variant identification using a genome tiling array and genotyping probes

Abstract

With over 5.5 million deaths worldwide attributed to the respiratory disease COVID-19 caused by the novel coronavirus SARS-CoV-2, it is essential that continued efforts be made to track the evolution and spread of the virus globally. The authors previously presented a rapid and cost-effective method to sequence the entire SARS-CoV-2 genome with 95% coverage and 99.9% accuracy. This method is advantageous for identifying and tracking variants in the SARS-CoV-2 genome compared with traditional short-read sequencing methods which can be time-consuming and costly. Herein, the addition of genotyping probes to a DNA chip that targets known SARS-CoV-2 variants is presented. The incorporation of genotyping probe sets along with the advent of a moving average filter improved the sequencing coverage and accuracy of the SARS-CoV-2 genome.

Keywords: COVID-19; SARS-CoV-2; bioinformatics; genotyping; pandemic; screening; sequencing; tiling-array; viral genome.

Figure 1

Overview of tiling array protocol and bioinformatics analysis strategy; a. Overview of sample preparation and handling. Samples are converted fron RNA to cDNA using ligase, followed by amplification biotin-11-dUTP, followed by fragmentation, hybridization to the tiling array, staining with Cy3 streptavidin, and imaging using a custom confocal scanner; b. Overview of base calling algorithm. Raw intensity measurements are ranked in decreasing order for each position probed. Values are used to calculate the difference (d) and differential (D) for each position. We make a ML model using these parameters and assign base calls. Resequenced data is aligned with the reference genome and used for further analysis; c. replacement iteration 1 looks at all variant calls. The variant calls made by the genotyping probeset may or may not match the call made by the genome tiling array. The Q scores of the two reads are compared. If the genotyping probeset confirms the variant call (ie genome tiling array = genotyping array) the call confirms the variant and the Q-score is replaced. If the variant calls between the two reads differ, the call and read with the higher Q-score is used to make the base call. All known variants at the time of array fabrication are accounted for; variants called by the genotyping array are automatically at ‘site of known variant’. If any calls are non-calls (ie Q-score < 20), when the genotyping probset Q-score is >20, calls are automatically replaced; d. replacement iteration 2 looks at all calls after replacement iteration 1. When he genotyping probeset makes higher quality base calls than the genome tiling array equivalent, the Q-score is replaced. The genotyping probeset can make correct reference or variant calls. If the Q-score is higher from the genotyping probeset and it confirms the reference or variant, base call and Q-score from the genotyping probeset is used. This raises the overall quality and confidence of the calls. ‘Confirmation’ refers to whether the genotyping probeset call matches the reference or the variant call made by the genome tiling array. ‘Non-confirmed’ indicates calls that do not match reference or known variant call; e. The MA-filter is constructed using the final Q-scores after replacement iteration 1 and 2. For each base call, the average Q-score of the 12 reads before and after are averaged to determine the MAQ at the position. The average of all MAQ – 2x standard deviation is used as the cutoff. All calls with a Q-score below 30 are entered into the filter. Of the intermediate quality calls, calls with a MAQ below MAQ – 2x standard deviation is designated non-call and the base call is not used.

Figure 2

Representative schematic of tiling array design and probe design for the genome tiling array and genotyping probe set; a. A small section of the genome tiling array. Each square is a feature covered with a specific 25mer oligonucleotide probe for a specific position for target genome. A probe set is a grouping of four features, one for each potential base. Highlighted in red are the probe sets for position 28122 and 28133. Probe sets are roughly organized in the order of base pair position; b. conceptual representation of tiling array probes for positions 28132, 28133, 28134. Each feature consists of a single 25mer probe with one of 4 possible bases at the 13^th position indicated by the green box which is the resequenced position. The base highlighted in blue is the reference base call. Others indicated in red are potential variants. Flanking 12bp on either side exactly match the reference base calls; c. Representative section of genotyping probe set design. The layout design, organization of probe sets, and 25mer design concept are identical to genome tiling array. Highlighted in red are features within the probe set for position 241, a site of known variant; d. conceptual representation depicting difference between probe sets for reference 241 (top) and variant 241* (bottom). Bases highlighted by the green box at the 13^th position indicate the probed position. The base highlighted in blue is the reference base call. For position 241, the reference base call is C while in 241* the reference base call is T. All other positions of the 25mer probes exactly match the reference base calls.

Figure 3.. Scatter plot, WY64 Q score and MAQ with variant call breakdown after MA-filter.

A scatter plot generated in R using ggplot2 displaying all reads made by the full genome tiling array reads with incorporation of genotyping probe-set data and MA-filter on sample WY64 between position 26 and 29834. Every call was assessed with a combination of base call and Q score. High quality reads with Q > Qth (20) are used to make calls. Dark gray circles represent base positions where Q > Qth and make a ‘Reference Call’. All calculated MAQ are overlaid as light gray. Brown dots are reads that have a low Q score where Q < Qth, are categorically ‘Non-calls’, and excluded to make base-calls. Q scores of ‘Variant Calls’ are identified as larger red circles, where the final base call made by DNA chip after replacement and MA-filter is not reference and the Q > Qth. Within ‘Variant Calls’, a blue overlap indicates calls that are removed by the MA-filter, which takes all variant calls with a Q score between 20 and 30 and removes any with a MAQ lower than the MA-threshold. Within ‘Variant Calls’, a green overlap indicates true variant calls identified and verified with short-read sequencing data. Additional colorized scatter plots for the remaining samples can be found under Supplemental Figures S1–S8.

Figure 4.. Venn Diagram, WY64 Variant Calls after MA-filter.

A Venn diagram depicting the categorical breakdown of calls made by the DNA chip on sample WY64. Variant Calls made by the DNA chip are contained within the red circle. Any circle or number outside of red implies that the call is a reference call. True variants are color coded in green and are confirmed by short-read sequencing data. The blue circle indicates variant calls that are filtered out by the MA-filter, converted to non-calls, and omitted from making base calls. Overlap between green and red circles indicates true variants correctly called by the DNA chip and verified through short-read sequencing data. Overlap between green and blue circles indicate an improperly removed true variant by the MA-filter where the read is situated in a local region of low MAs indicated by a low MAQ below the MAQ threshold. Ideally both green and blue circles should be contained within the red circle without overlapping each other, implying that all true variants and reads filtered out by the MA-filter were correctly identified as variant calls. Additional Venn diagrams for the remaining samples can be found under Supplemental Figures S9–S16.