key: cord-0907782-jmwdsf18
authors: Alizon, S.; Selinger, C.; Sofonea, M. T.; Haim-Boukobza, S.; Giannoli, J.-M.; Ninove, L.; Pillet, S.; Vincent, T.; de Rougemont, A.; Tumiotto, C.; Solis, M.; Stephan, R.; Bressollette-Bodin, C.; Salmona, M.; L'Honneur, A.-S.; Behillil, S.; Lefeuvre, C.; Dina, J.; Hantz, S.; Hartard, C.; Veyer, D.; Delagreverie, H. M.; Fourati, S.; Visseaux, B.; Henquell, C.; Lina, B.; Foulougne, V.; Burrel, S.
title: Epidemiological and clinical insights from SARS-CoV-2 RT-PCR cycle amplification values
date: 2021-03-17
journal: nan
DOI: 10.1101/2021.03.15.21253653
sha: a421ce37ab5a2b7296ae369341c96da3e7285d6e
doc_id: 907782
cord_uid: jmwdsf18

The SARS-CoV-2 pandemic has led to an unprecedented daily use of molecular RT-PCR tests. These tests are interpreted qualitatively for diagnosis, and the relevance of the test result intensity, i.e. the number of amplification cycles (Ct), is debated because of strong potential biases. We analyze a national database of tests performed on more than 2 million individuals between January and November 2020. Although we find Ct values to vary depending on the testing laboratory or the assay used, we detect strong significant trends with patient age, number of days after symptoms onset, or the state of the epidemic (the temporal reproduction number) at the time of the test. These results suggest that Ct values can be used to improve short-term predictions for epidemic surveillance.

Molecular testing is a key component of screening policies to control the spread of infectious diseases and the SARS-CoV-2 pandemic has led to an unprecedented testing rate using reverse transcription polymerase chain reaction (RT-PCR) assays. (1) . In clinical and public health practices, RT-PCR results are qualitative for viral respiratory disease diagnostics, with reports such as 'positive', 'negative', 'uninterpretable', and, sometimes, 'weakly positive'. These are based on the cycles threshold, also referred to as crossing point or crossing threshold (here denoted C t ), which corresponds to the number of PCR amplification cycles required for the fluorescent signal to rise above a positive threshold. In theory, the more abundant the genetic target in the sample, the fewer the amplification cycles required to detect it. This is why numerous studies on SARS-CoV-2 rely on C t values to assess transmissibility (2) or study infection kinetics (3) . Here, we present a cross-sectional analysis of SARS-CoV-2 RT-PCR tests performed on 2,220,212 individuals in France between January 21, 2020, and November 30, 2020 ( Figure 1 ). Few studies analyze C t values at a population level. One explanation for this matter of affairs is that these are known to suffer from several, potentially strong, biases. First, sample type and sampling quality directly affect the amount of genetic material available. Second, the RT-PCR assay matters. Even the quality of the reagents used may have a significant effect on the number of amplification cycles required to achieve the same level of fluorescence for the same amount of target genetic material. Combining data from different laboratories helps to control for these sources of variation in statistical analyses. Furthermore, the larger the dataset, the more we can detect small statistical trends even after having controlled for non-informative variables. In our analysis, we studied tests from individuals aged between 1 to 89 years old. We did not take into account tests for which key variables such as patient age, patient sex, laboratory geographical department, qualitative result, or RT-PCR assay used were unknown. Note that one test could provide more than one C t value if containing probes targeting multiple viral genes. According to the national guidelines (4), it is recommended to focus on the most sensitive target to categorize levels of viral excretion. After removing the biologically unrealistic C t values that were lower than 10 or larger than 45, the 95% confidence interval (95CI) of the remaining values was [16.89; 35 .56] (Table S1 ). These correspond to 793, 479 tests from the same number of individuals. We used a linear regression model to explore how C t values can be explained by the following variables: patient age and sex, the number of days since the onset of the symptoms (if known), the clinical sampling site (if known), the sampling facility (if known), the RT-PCR assay used, the target gene, the test's qualitative result, the sampling date, the temporal reproduction number of the epidemic at the sampling date (denoted R t and estimated using the EpiEstim method (5)), and a control variable. The latter corresponds to the last digit of the patient anonymity number and is expected to be independent of the C t value. We also included in the model an interaction term between sampling date and R t . For this analysis, we excluded C t values from internal controls. Univariate analyses are extremely sensitive to heterogeneity in the dataset. For instance, the age distribution from patients sampled in aged care homes is different than that from city screening facilities, and analyzing the 'sampling facility' factor alone could yield misleading result. This is why the analysis used here is multivariate and considers all the factors listed above simultaneously. Overall, the linear model explained 38.8% of the variance in C t values, and the residuals were normally distributed ( Figure S2A ). All the factors except the temporal reproduction number were significantly associated with C t values using a classical 5% p-value criterion in an analysis of variance (ANOVA) with type II sums of squares. Even for the control variable, the p-value was 0.013 and patients with final digits 1 and 3 had C t values slightly lower (-0.19 and -0.17 cycles) than patients with a 0 final digit. We therefore set our significance thresholds to 5% of that of the control variables, i.e. 6.5 × 10 −4 to analyze main effects (Table 1) . The intercept of the linear model indicates the average C t value for a positive test performed with the reference assay, and all the other factors being set to their reference value. Its magnitude (19.1 cycles) is in line with clinical practice. The importance of the noise in the dataset is illustrated by the strong effect of the testing laboratory, as well as the RT-PCR assay used (Supplementary Figure S1) . Despite this high level of noise, we detected a strong effect of the qualitative result (Figure 2A ), with C t differences that were even larger than that from the laboratory effect. We also found a slightly significant difference of -1.81 cycles between the most common type of samples (nasopharyngeal) and that performed on other clinical sampling sites (mostly lower respiratory tracts, but also feces or saliva). This is likely because the latter tests were performed in patients with more severe symptoms. The effects associated with the number of days since symptoms onset was particularly strong. For the 8.5% of the participants for whom the number of days between symptoms onset and testing dates was known, we found that the C t gradually increases over the reported range with a maximum difference of 5.73 cycles ( Figure 2B ). The effect of sex had the same order of magnitude as that of the control variable and could, therefore, be treated as non-significant. Conversely, the age factor had a strong effect with a decrease of 0.541 cycles per year ( Figure 2C ). The target gene of the RT-PCR assay used also yielded a slightly significant effect. The C t values obtained when using a probe targeting the ORF1 and S regions of the virus genome were significantly higher than when using the N gene, which was the genomic region of reference in the model ( Figure  2D ). This effect is consistent with the life-cycle of the virus. As stressed by (6) , since coronaviruses are (+)ssRNA viruses, they use the same RNA matrix for replication and transcription, both being amplified by diagnostic assays. Furthermore, Coronaviridae transcripts can produce subgenomic mRNAs that lack part of the genome (7) . As a consequence, and as shown in cell cultures (8) , genes on the 5' end of the genome are under-represented. This is consistent with our result where assays targeting the gene on the 3' end (the N gene) tend to have lower C t than assays targeting genes on the 5' end (the ORF1 and S genes). Note that an alternative explanation could be that some probes target more conserved areas of the SARS-CoV-2 than others (9) . Finally, we found that C t values decreased with time (-0.797 cycle per day), but this effect was non-linear ( Figure 2E ). This could be due to the strong variation in testing efforts in France ( Figure 1A ), but also to variations in the epidemic trend. Indeed, although the R t (inferred from hospitalization data using the EpiEstim method (5)) was not found to be significant, the interaction between the sampling date and R t was nearly significant ( Figure 2F ), suggesting that a temporal analysis could yield additional insights. The existence of a correlation between the C t values of the tests performed in a population and R t is consistent with population dynamics theory, which predicts that in an expanding population of infected individuals, the 'age' of the infections, i.e. the number of days post-infection, is skewed towards lower values (10) . Since C t values have been reported to increase over the course of an infection (3), which we confirm with this analysis ( Figure 2D ), it has suggested that these values could be used as an early signal to predict R t (11). To investigate this question, we focused on screening data collected in the general population from individuals aged from 5 to 79. We estimated the median and skewness values of the daily distribution of the C t values. To correct for potential confounding factors, these were adjusted using a linear model (see the Supplementary Methods). We analyzed the temporal correlation between the time series with a 7-day rolling average of this median, skew, and R t (Figure 3 ). For the median C t value, we found a significant correlation with R t that was maximized for a 6 to 7 days delay ( Figure S2 ). This is consistent with R t being calculated using data from ICU-admissions, which occur with a median of 14 days after infection (12, 13) , and RT-PCR screening data being obtained earlier in the infection. To further assess the usefulness of C t data, we used ARIMA models to predict R t dynamics over 7 days. We compared models without any exogenous data, to models that also included exogenous time series (either median or skewness of estimated C t values distribution, or the fraction of positive tests (1)). As expected, the prediction error made using only endogenous data (R t ) was low in periods where R t variations were limited. Furthermore, we found that adding exogenous data improved the prediction, especially when strong shifts in R t were occurring ( Figure 3 ). C t values (green and cyan dots) tended to provide a better reduction in the error of the prediction than the ratio of positive tests (purple). This analysis of a large national database of RT-PCR tests performed in the context of a major epidemic confirms that population-level C t values are noisy since even a linear model that features 91 degrees of freedom does not explain the majority of the variance. However, owing to the law of large numbers, we detect several effects that are in line with biological observations and with virological properties. For instance, our finding that C t values decrease as a function of the number of days after symptoms onset is consistent with longitudinal follow-ups (3) . Similarly, the difference we detect between the virus gene targeted by the RT-PCR assay used can be interpreted in the light of known differences in mRNA copy numbers between genes depending on their distance from the 3' end (8) . Concerning the link between age and C t values, although there are some mechanistic interpretations as to why virus load would increase with age (14), the evidence was mixed, with some studies reporting a decreasing trend (15) and others not (16, 17) . Here, using a multivariate approach on a large dataset allows us to unravel a strong and significant decrease of C t values with age. A limitation our study is that although our dataset stands out by its size and its level of details, it is restricted to a single country where testing effort varies, both on a temporal and on a spatial scale ( Figure  1 ). Performing similar analyses for other populations can therefore be particularly informative. A promising output of this analysis is the possibility to use C t values as an early signal to detect changes in epidemic behavior, e.g. R t values. Indeed, our most robust descriptors of the epidemics originate from hospital-admission data, but these have a significant delay with the status of the epidemic since patients are hospitalized 2 weeks after infection (12, 13) . The ratio of positive tests performed in the population of interest can, in theory, provide earlier insights but it suffers from strong sampling biases. We show that accounting for population-based C t values variations can improve R t predictions on a 7-days period. This corroborates earlier hypotheses (11) but also reveals the importance of analyzing a large-enough dataset to filter out the important amount of noise in these values. Overall, these results call for better integration of C t values in national public health policies to monitor epidemics caused by SARS-CoV-2 but also other human viruses, especially since these data raise less ethical concerns than other sources of data such as mobility data. 5 All rights reserved. No reuse allowed without permission.

(which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted March 17, 2021. 8 All rights reserved. No reuse allowed without permission.

(which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted March 17, 2021. All rights reserved. No reuse allowed without permission.

(which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. , the median (in green) and skewness (in cyan) of the daily C t residual distribution, and R t (in black). The bottom panels show the error made by a prediction using only R t data (red dots) and the potential improvement made by including exogeneous data.

10 All rights reserved. No reuse allowed without permission.

(which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. 

This study was approved by the Internal Review Board of the Centre Hospitalier Universitaire de Montpellier. It is registered at ClinicalTrials.gov under the identifier NCT04738331.

The final data set analyzed along with the R scripts will be made available upon publication. (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. Table S1 . For each individual, we retained only the earliest sample therefore analyzing 1, 129, 437 C t values.

13 All rights reserved. No reuse allowed without permission.

(which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted March 17, 2021. ; Table S1 : Description of the dataset variables. For real numbers, we show the median and the 95% confidence interval. For categorical variables, we either indicate the number of factors or the number of occurrences n for each factor. 

The temporal reproduction number (denoted R t ) was computed on the COVID-19 hospital admission incidence time series established by the national public health agency (Santé Publique France) and accessible at this website. Because of strong daily variations (especially on week-ends), we first transformed the time series using a 7-days rolling average. We then used the EpiEstim method (5) and the eponym package in R (18). Earlier studies have reported that, for patients who develop severe symptoms, the median time between infection and hospital admission is in the order of 14 days (12, 13).

All the analyses were performed in R version 3.6.3.

The linear model was performed in individuals from 1 to 90 years old. We also removed C t values associated to internal controls because they caused the distribution of the residuals to be non-Gaussian. Quantitative factors, namely R(t), date, and age were scaled and centered. The model was formulated as follows in R: modele_general = lm(Ct~Rt * date + age + sexe + target_gene + assay_PCR + id_lab + symptoms + result + sample_type + sampling_facility + control_variable, data)

The adjusted R 2 of the model was 38.8% and the distribution of the residuals Gaussian ( Figure S2A ). The model was analyzed using an ANOVA assuming type-II errors because of the unbalanced nature of the data set. All the variable were found to be extremely significant (p-value < 10 −6 ), except for R t (p-value of 0.68) the control variable (p-value of 0.0131), sampling facility (p-value of 0.0135), sex (p-value of 0.0021), interaction between R t and date (p-value of 0.000842).

To analyze the time series of reproduction number and C t values we restricted the data to tests performed after July 1, 2021, in a screening context, using nasopharyngeal swabs, in individuals aged from 6 to 80. These assumptions were made such that C t values would reflect the state of the ongoing epidemic. Finally, we excluded values from internal control genes. Overall, we analysed 234,782 C t values from 110,227 individuals. To correct for other potential confounding factors, we first performed a linear model. The analyses were performed on the residuals of the linear model. For each day, we computed the median and the skewness of the distribution. Skewness was computed using the formula n (n−1)(n−2) i

, with x i the individual values and n the total number of values. We then computed a 7-day rolling mean to buffer daily variations. The temporal reproduction number was computed as indicated above. We also analyzed a 7-day rolling mean. Cross-correlation function analyses were performed using the ccf function in R.

We used ARIMA models (implemented in the R package forecast) to predict the hospital-admission temporal reproduction number (R t ) from past observations. For each date, predictions were evaluated in terms of the mean absolute percentage error (MAPE) for horizons of 7 days in the future based on coefficients learned from past data starting on July 29th, 2020. More precisely, for each date we compared the temporal reproduction number data {D k ; k = 1, . . . , 7} with the seven-day model forecast {F k ; k = 1, . . . , 7} by

We considered 4 types of data:

1. past R t data (i.e. endogenous data), 2. quartiles from C t residuals (to remove biases), 3 . skewness from Ct residuals, 4 . national positive test ratio collected from https://covid.ourworldindata.org/data/.

The residuals of C t were obtained from the following linear model:

lm(Ct~age + target_gene * assay_PCR + id_lab, data)

Our goal was to see if adding exogenous data, i.e. C t values and proportion of positive tests, increased prediction precision. Models were tuned with the auto.arima function, and untuned models were run with p = 9, d = 2, and q = 0 as default parameters, based on the cross-correlation analysis between R t and C t time series ( Figure S2 ). Prediction improvement by models using in addition exogenous data relative to model using past values of R t only was defined by ∆MAPE = MAPE Rt − MAPE Rt+exo

To visualize the effect of the different factors, we first perform a linear model to correct for the effect of confounding factors.

We used the general model described above but removed the effect 'control_variable', as well as the effects 'target_gene', 'sample_type', 'symptoms' and 'sampling_facility' which had little effect and were sometimes lacking for many participants. For each figure, we also removed the main factor of interest (depending on the panel). The C t value itself was generated using the predict.lm function in R.

All rights reserved. No reuse allowed without permission.

(which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. Figure S1 : Effect of the RT-PCR assay used on the estimated C t value as a function of the targeted genomic area. Only tests with at least 1,000 C t values are shown.

All rights reserved. No reuse allowed without permission.

(which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. Figure S2 : Distribution of the C t residuals value. A) For the linear model used for the main analysis shown in Table 1 , B) for the linear model used to generate residuals for the R t time series analysis. q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q −0. q q q q q q q q q q q q q q q q q q q q q q q q q q Figure S3 : Cross correlation functions between R t and A) the median or B) the skewness of the C t residuals distribution. The blue shaded areas show the non-significant values (with a 95% threshold) and the red vertical dotted lines the lag with the highest significant correlation. Note that the lag is smaller for the skewness than for the median of the distribution.

18 All rights reserved. No reuse allowed without permission.

(which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. Figure S4 : Predicting temporal reproduction number (R t ) from time time series with untuned ARIMA parameters p = 9, d = 2, q = 0 The top four panels show the the 7-days rolling averages of the time series of the ratio of positive tests (in purple), the median (in green) and skewness (in cyan) of the daily C t residuals distribution, and R t (in black). The bottom panels show the error made by a prediction using only R t data (red dots) and the potential improvement made by including exogeneous data.

All rights reserved. No reuse allowed without permission.

(which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. 

Avis du 25 septembre 2020 de la Société Française de Microbiologie (SFM) relatif à l'interprétation de la valeur de Ct (estimation de la charge virale) obtenue en cas de RT-PCR SARS-CoV-2 positive sur les prélèvements cliniques réalisés à des fins diagnostiques ou de dépistage

Matrix population models: construction, analysis and interpretation