Using intrahost single nucleotide variant data to predict SARS-CoV-2 detection cycle threshold values

Abstract

Over the last four years, each successive wave of the COVID-19 pandemic has been caused by variants with mutations that improve the transmissibility of the virus. Despite this, we still lack tools for predicting clinically important features of the virus. In this study, we show that it is possible to predict the PCR cycle threshold (Ct) values from clinical detection assays using sequence data. Ct values often correspond with patient viral load and the epidemiological trajectory of the pandemic. Using a collection of 36,335 high quality genomes, we built models from SARS-CoV-2 intrahost single nucleotide variant (iSNV) data, computing XGBoost models from the frequencies of A, T, G, C, insertions, and deletions at each position relative to the Wuhan-Hu-1 reference genome. Our best model had an R2 of 0.604 [0.593-0.616, 95% confidence interval] and a Root Mean Square Error (RMSE) of 5.247 [5.156-5.337], demonstrating modest predictive power. Overall, we show that the results are stable relative to an external holdout set of genomes selected from SRA and are robust to patient status and the detection instruments that were used. This study highlights the importance of developing modeling strategies that can be applied to publicly available genome sequence data for use in disease prevention and control.

Fig 1. Histograms showing the distributions of the 20 most common variants in the dataset.

Fig 2. Histograms showing the distributions of Ct values in the dataset.

Ct values were sorted into bins of 3 with an inclusive lower bound and exclusive upper bound.

Fig 3. Schematic of the matrix generation for the XGBoost models.

Pileup files were generated by aligning reads against the reference genome, and the iSNV frequencies of A,T,G,C, insertions, and deletions at each position were computed per genome and normalized with respect to read depth. Normalized iSNV values and the one-hot-encoded clinical detection instruments were used to create the matrix that was used to generate the XGBoost models.

Fig 4. Scatterplots of predicted versus actual Ct values.

Scatterplots were constructed for models trained and tested on the following instruments: a) All instruments, b) Alinity, c) Panther, d) Cepheid. Points are colored by variant with samples of the 10 most frequently occurring variants colored via the key shown in the right, and samples of other variants colored gray. The line y = x is shown across the center diagonal of the figure for reference. Data are from a single fold.

Fig 5. Confusion matricies for models by instrument using a single fold.

Model predictions were binned into Ct value ranges of 3 cycles with an inclusive lower bound and exclusive upper bound. Coloring and values in each cell represent the fraction of the actual Ct values predicted in the given interval. Empty cells with no predictions or actual values in that range are gray. Confusion matrices were constructed for models trained and tested on the following instruments: a) All instruments, b) Alinity-only, c) Panther-only, and d) Cepheid-only.

Fig 6

Dot plot depicting the average Ct value (A) and XGBoost feature information gain (B) for each position and character used by the all-instrument model. Each base at a given position is colored according to the key. Genomic positions correspond to the SARS-CoV-2 Wuhan-Hu-1 reference genome. For each position, only genomes where ≥40% of the characters in the column corresponded to a given nucleotide were used to generate the average gain in order to reduce noise in the image. Additionally, only statistically significant bases are included, significance was computed based on the 95% confidence interval of the average Ct value of genomes with a given base and those without. No INDEL features met this significance requirement. The spike protein corresponds to genomic coordinates 21563–25384.

Fig 7

The A) R² (red line), and B) RMSE (blue line) for a holdout set of Omicron genomes using models trained on increasing percentages of Omicron genomes in the training set. The green dashed line depicts the R² and RMSE for the training set, which contains no Omicron genomes.