Genomic epidemiology offers high resolution estimates of serial intervals for COVID-19

Simulation study

In the simulation study, with results presented in Fig. 2, we simulate an influenza-like outbreak using the outbreaker package in R²¹, that simulates transmission trees with associated symptom onset times and pathogen genomic sequences. To explore the impact of incomplete case sampling, we run experiments in which we mask an increasing proportion of infected hosts from the serial interval estimation procedure. That is, we run our method with 100%, 90%, …, 10% of infected cases only, chosen at random. The outbreaker simulator assumes a simple SIR outbreak in which symptom onset and onset of infectiousness occur together, and hence the serial interval is exactly equal to the generation time. We simulate the outbreak with a susceptible population size of 1000, R₀ = 2, mutation rate (per site per day) 10⁻⁴ and case importation rate 0.01 per day. The true serial interval distribution used for simulation is Γ(μ = 4.5, σ = 2). This results in an outbreak of final size 807.

In the simulation study, we provide our method with the symptom onset time and pathogen genomic sequence from each sampled case, but hide the true transmission tree. For each experiment, in which a proportion p = 1.0, 0.9, . . ., 0.1 of infected cases are included, we first construct a transmission cloud of feasible transmission pairs, by taking all pairs with genomic distance <0.0001 (see ‘Sampling transmission networks’ below). This distance is chosen to be approximately as stringent as our Victorian COVID-19 analysis, given the simulated mutation rate. We sample a larger number of transmission networks (1000) than in the Victorian analysis (100), to reflect the larger population size, that may lead to a larger space of plausible networks. From the sampled networks, we obtain MAP estimates of the model parameters (μ, σ, π, w), using the same approach as for the Victorian COVID-19 outbreaks, described in full below. We use beta distributed prior distributions on π and w that represent moderate knowledge of the sampling rate, as this is the situation in which our method is intended to be most suitable. That is, the prior mean is set between 90–110% of included proportion p (varied across experiments), with the prior standard deviation set to 0.1. We note that p is not exactly equivalent to either π or w, the relationship will depend on the structure of the sampled network. We do not use priors for μ or σ.

Victorian SARS-CoV-2 whole-genome sequencing and clustering

The remainder of the Methods describes the analysis of the Victorian SARS-CoV-2 clusters. All samples were routinely sequenced as part of public health operations. We provide a brief description of the sequencing and genomic clustering procedure below, for full details, see ref. ⁶. Whole RNA was extracted from samples obtained from nasopharyngeal swabs, and positive samples for SARS-CoV-2 were identified by RT-qPCR. Positive samples underwent tiled amplicon sequencing using ARTIC primers (v1 and v3)^36,37. Amplicons were prepared for sequencing on an Illumina sequencer using the NexteraXT library prep protocol following the manufacturer’s instructions. Sequencing reads were mapped to the Wuhan-Hu-1 reference sequence (Genbank MN908947.3 [https://www.ncbi.nlm.nih.gov/nuccore/MN908947]) and consensus sequences generated using the iVar pipeline³⁸. Consensus sequences were kept for downstream analyses if they met the following criteria: ≥95% genome recovered, ≤25 SNPs from the reference genome, and ≤300 ambiguous bases. The data comprise genomic sequence, sequence sampling date and symptom onset date for each sampled case, in addition to some epidemiological information for wave 2 as described below.

Wave 1 clustering

A total 1242 samples from 1075 patients were sequenced during the wave 1 study period: 6 January–14 April 2020. The wave 1 genomic clustering procedure is as described in ref. ⁶. Of 903/1242 samples passing initial quality control and de-duplication, 737 were identified as belonging to a genomic cluster, for a total of 76 clusters with at least two cases. The clustering of the samples was based on a maximum-likelihood phylogenetic tree containing a single sequence per patient built using IQtree³⁹ applying a GTR+Γ4 substitution model. Using the ClusterPicker tool⁴⁰, clusters were defined as having at least two samples with the inferred ancestral node having at least 95% bootstrap support and the maximum distance within the cluster of 0.0004 expected substitutions/site. The 76 clusters had strong concordance to epidemiologically-defined clusters (a median of 100% of cases in each epidemiologically-linked group were in a single genomic cluster, but genomic clusters were wider: with a median of 43% of cluster cases in a single epidemiological group)⁶. We focus on all genomic clusters with at least 15 cases for the serial interval analysis. After removing sequences with <95% genome recovered or with missing symptom onset time (three cases), the wave 1 data comprise a total of ten wave 1 clusters containing a total of 312 cases.

Wave 2 clustering

Even when pathogen samples are routinely collected, genomic clustering can have computational challenges and lead to large, uncertain clusters, particularly when the transmission is widespread in a community and many cases are sampled. In this context, the concept of what constitutes a cluster is less clearly defined, and tree- or sequence-based clustering can lead to infeasibly large clusters. This is especially true for SARS-CoV-2, for which sampled genetic diversity has been described as remarkably low⁴¹, even more so during periods of lockdown and limited international travel⁴², leading to challenges in phylogeny building and sub-classification⁴³. Public health experts may also be interested in transmission not simply among genetically linked isolates, but within particular locations or groups of people, for example, schools and hospitals. Driven by these factors, for the wave 2 analysis, we derive clusters epidemiologically rather than genomically: choosing all cases associated with particular exposure sites defined by public health.

A total 15,665 samples from 14,075 patients were sequenced during the wave 2 study period: 1 June–28 October 2020. Epidemiological data, including detailed demographic, risk factor and contact tracing data, were collected for each case through an interview conducted by the Victorian Department of Health, and these data were used to determine locations of possible transmission, known as exposure sites. In our primary analyses, we use only the exposure site information (which cases attended which sites) to define the wave 2 clusters, we do not use any demographic or contact data. We perform a secondary analysis in which we re-estimate the serial intervals without the viral sequences but given full knowledge of the contact data.

Of the 10,642/15,665 samples passing quality control and de-duplication, 7116 were identified as associated with at least one exposure site. Of these samples, we remove 1371 cases with missing symptom onset time. This provides a total of 623 clusters with at least two cases. Recall that cases may be associated with more than one exposure site, and so may be included in more than one cluster. In our primary analysis, we focus on ten exposure site clusters—we select the largest two clusters associated with each of the aged care facilities, healthcare facilities, housing, schools and meat packing/meat processing plants. We perform additional analysis on all wave 2 clusters with at least 15 cases. We remove a further 2 clusters which contain more than 15 cases but have a limited signal of transmission (five or fewer plausible infector-infectee pairs). This results in a total of 94 wave 2 clusters comprised of a total of 3875 cases.

Phylogenetic reconstruction

Phylogenies are produced from the entire set of sequences in each wave, using a custom workflow for building fast SARS-CoV-2 trees available at github.com/MDU-PHL/kovid-trees-nf. Briefly, sequences are cleaned to remove sites with >5% missing calls and de-duplicated with GOALIGN⁴⁴. An approximate maximum-likelihood tree is built using FastTree⁴⁵ and the branch lengths are optimised with RAxML-NG⁴⁶. Finally, GOTREE⁴⁴ is used to repopulate the tree with duplicate sequences.

Sampling transmission networks

In order to estimate cluster-specific serial intervals from viral sequences and symptom onset times in a cluster (i.e. in the absence of contact tracing data), we require an approach for obtaining putative transmission pairs. Given that there will usually be considerable uncertainty in who infected whom, we sample a set of plausible transmission networks, accounting for indirect transmission and perform parameter estimation across the entire set.

We build a pairwise genetic distance matrix between all aligned sequences in a cluster using the ape package in R and the TN93 model of evolution (this was found to generate very similar distances to the GTR+Γ4 model used to generate the wave 1 clusters, see Supplemental Materials Fig. S1). We then create a ‘transmission cloud’, that is, a set of all plausible transmission pairs in the cluster. We define a plausible (possibly indirect) transmission pair as two cases in the same cluster who:

1. 1.

   have observed interval between symptom onsets 0 < T ≤ 35 days

2. 2.

   have pairwise genomic distance G < 1.1/29903 (where 29,903 is the alignment length).

We apply a much more stringent genomic distance criterion (corresponding to ~2 generations of infection) than the onset distance criterion (corresponding to ~7 generations), to allow the genomic data to be relatively more informative than the time-based data, and to minimise input bias towards short serial intervals. Note that our approach may result in a case having multiple plausible infectors, or no plausible infectors: in the second scenario, the case would be considered as an importation to the cluster from an unsampled case.

From the transmission cloud, we sample a plausible transmission network by sampling an infector for each infectee from among their list of plausible infectors. For infectee j, a plausible infector i is selected with higher probability if the genomic distance and/or the difference in symptom onset time between i and j is lower. The sample weighting s(i, j) is given by

where G_i,j is the genomic distance between i and j, and T_i,j is the difference in symptom onset time between i and j. \({G}_{\max }\) and \({T}_{\max }\) are the maximum genomic and time distance, respectively, among all plausible pairs in the cluster; their inclusion causes distances G_i,j and T_i,j close to zero to result in higher sampling probability s_i,j. The weighting is normalised to sum to 1 for each infectee j. We repeat this process N = 100 times to obtain a set of 100 plausible transmission networks per cluster.

As we do not assume that all cases have a sampled ancestor within their cluster, a sampled transmission network is not required to be comprised of a single transmission tree; there can be several distinct sub-trees, each spawned by an unsampled case. The result is that we only perform serial interval estimation on those inferred transmission pairs with sufficient confidence.

Serial interval estimation

After sampling a set of plausible transmission networks per cluster, we proceed to estimate the serial interval in each network, and finally combine these for a single cluster estimate. Our approach seeks to estimate the parameters of the true underlying serial interval distribution, that is, the serial interval arising from direct transmission between a pair of individuals. However, we account for the fact that observed serial intervals in data with incomplete sampling may comprise a mixture of direct transmission (i → j), indirect transmission (i → x → j, for any number of unsampled x) and coprimary transmission (x → i, x → j, for unsampled x).

In this description, we assume that a plausible transmission network has already been sampled from the data as above, and so transmission between a sampled infector-infectee pair (i, j) is considered certain (even if that transmission may be direct, indirect or coprimary). We apply the methodology described in this section independently for each sampled network, in order to incorporate uncertainty in the network.

Let \({T}^{{\prime} }\) be the true serial interval distribution of the disease under consideration, so that \({T}_{ij}^{{\prime} }\) is the serial interval arising from direct transmission between any case i and case j. We assume that \({T}_{ij}^{{\prime} } \sim \Gamma (\mu,\sigma )\), for mean μ and standard deviation σ. Note that we assume, therefore, that the serial interval is strictly non-negative. This is the quantity we are interested to make inference about, by estimation of μ and σ.

Now, let T_ij denote the observed time interval between the symptom onsets of a particular case i and case j. If there is direct transmission from i to j then \({T}_{ij}={T}_{ij}^{{\prime} }\). However, if we allow for under-reporting in the data, i.e. there may have been unsampled intermediate cases between i and j, then T_ij is the convolution of multiple serial intervals. It is also possible that both i and j were infected by the same unsampled host, which we call coprimary transmission. Although we allow for a theoretically infinite number of unsampled intermediate cases between cases i and j, in the coprimary case, we assume that a single unsampled host must have infected both i and j.

More formally, we assume that the serial interval distribution is a mixture of two transmission paths: coprimary and non-coprimary (direct or with unsampled intermediates).

Non-coprimary transmission

Sampled infector i and infectee j are separated by m unsampled intermediate hosts, 0 ≤ m < ∞ (m = 0 corresponds to direct transmission). Then, the observed interval T_ij is the sum of m + 1 i.i.d. gamma-distributed intervals, \({T}_{ij} \sim \Gamma ((m+1)\mu,\sqrt{(m+1)}\sigma )\). The value of m is unknown but can be considered to follow a Geometric distribution with success probability π, where success corresponds to sampling an infected individual and m is the number of failures. The parameter π then represents a pseudo-sampling probability, for hosts descending from a sampled infector only. Taking this into account, the marginal distribution of T_ij will be a Compound Geometric Gamma distribution, CGG(μ, σ, π). We assume that a proportion w of observed serials are of this type.

Coprimary transmission

In the coprimary case, sampled infector i and infectee j were in actuality both infected by the same unsampled case x. Since the serial intervals \({T}_{xi}^{{\prime} }\) and \({T}_{xj}^{{\prime} }\) are i.i.d. Γ(μ, σ), the observed interval T_ij is the strictly non-negative difference of two i.i.d. gamma distributions, \({T}_{ij} \sim \left|{T}_{xj}^{{\prime} }-{T}_{xi}^{{\prime} }\right|\). We refer to this as a Folded Gamma Difference distribution FGD(μ, σ). This is a folded version of the Gamma Difference distribution introduced in refs. ^47,48, with the simplification that both gamma distributions have the same parameters. The remaining proportion of 1 − w of observed serials are of coprimary type.

Likelihood

We wish to write down the likelihood of the model described above, given a dataset \({{{{{{{{\mathcal{D}}}}}}}}}_{c,\tau }={\{{T}_{{i}_{1},{j}_{1}},{T}_{{i}_{2},{j}_{2}},\ldots,{T}_{{i}_{n},{j}_{n}}\}}_{c,\tau }\) which contains the observed time intervals for all n sampled infector-infectee pairs in network τ for cluster c. Note that any given host can be an infector of multiple cases and an infectee as well, so some i and j may refer to the same individuals.

The log-likelihood of the model parameters given dataset \({{{{{{{{\mathcal{D}}}}}}}}}_{c,\tau }\) is given by:

where f_CGG and f_FGD are PDFs of the Compound Geometric Gamma and Folded Gamma Difference distributions, respectively.

Maximum a posteriori (MAP) estimation

Rather than maximising the likelihood expression directly to obtain estimates of the model parameters, we take a Bayesian approach to incorporate prior information on parameters π and w, and hence find maximum a posteriori (MAP) estimates. Given that there may be considerable uncertainty on the true transmission network as obtained from viral sequences and symptom onset times alone, but we often have good prior knowledge of the level of sampling in the population, the use of priors can help to avoid issues of identifiability in the model. We can imagine that, without the restriction of such priors, extremely short serial intervals with very low levels of sampling would provide a good fit to any data, even if we think this is impossible in practice.

We assume a beta distributed prior for both π and w, desirable as it is restricted to the range [0, 1]. The log prior distributions for π and w are given by:

Given these prior distributions and the log-likelihood expression introduced in Equation (2), the log posterior distribution is given by:

MAP estimates \({(\hat{\mu },\hat{\sigma },\hat{\pi },\hat{w})}_{c,\tau }\), conditional on a sampled transmission network τ for cluster c, are then found by maximising the log posterior density in Equation (5). In practice, we numerically optimise the log posterior density using the optim function in R.

What remains is to incorporate uncertainty in the transmission network τ. Rather than sampling a single transmission network τ from the transmission cloud of potential infector-infectee pairs for cluster c, we sample a set of networks τ₁, τ₂, …, τ_N. MAPs \({(\hat{\mu },\hat{\sigma },\hat{\pi },\hat{w})}_{c,{\tau }_{i}}\) are obtained for each network τ_i independently. When fewer than N MAP estimates are returned from the optimisation, as can occur when networks are randomly sampled which are not concordant with the assumed priors and so cause the numerical optimisation to fail, we sample additional networks until N MAPs are obtained.

The overall MAP for each parameter in cluster c is then calculated as the mean across all sampled transmission networks, for example:

and similarly for \({\hat{\sigma }}_{c},{\hat{\pi }}_{c}\) and \({\hat{w}}_{c}\).

Confidence Intervals

We obtain estimates of the standard error (\(\hat{\,{{\mbox{se}}}\,}\)) of MAPs \({(\hat{\mu },\hat{\sigma },\hat{\pi },\hat{w})}_{c,\tau }\), for a sampled network τ, using the inverse negative Hessian evaluated at the MAPs, as obtained from the numerical optimisation procedure. At confidence level α, this provides an approximate confidence interval for \({\hat{\mu }}_{c,{\tau }_{k}}\) of

and similarly for \({\hat{\sigma }}_{c,{\tau }_{k}},{\hat{\pi }}_{c,{\tau }_{k}}\) and \({\hat{w}}_{c,{\tau }_{k}}\).

In order to obtain confidence intervals for the overall cluster estimates \({\hat{\mu }}_{c},{\hat{\sigma }}_{c},{\hat{\pi }}_{c},{\hat{w}}_{c}\) i.e. over the space of all sampled transmission networks, we must take into account variation both within and between estimates. The variance of the estimator \({\hat{\mu }}_{c}\) (and equivalently \({\hat{\sigma }}_{c},{\hat{\pi }}_{c}\), and \({\hat{w}}_{c}\)) is derived with the law of total variance:

So, the first term incorporates the average uncertainty in each estimate of the MAP (for each sampled transmission network), and the second incorporates the estimate’s variability between sampled transmission networks.

This induces a confidence interval for the cluster-level estimates of

Inclusion and ethics

Data were collected in accordance with the Victorian Public Health and Wellbeing Act 2008. Ethical approval was received from the University of Melbourne Human Research Ethics Committee (study number 1954615.3).

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.