Data.

We analyzed the consensus sequences deposited in the GISAID database (4) with high quality and full lengths (number of bps ≈30,000). Four datasets are used for our investigation according to the collection date in GISAID database. The dates are 1 April 2020, 8 April 2020, 2 May 2020, and 8 August 2020. The list of GISAID sequences used is given in Dataset S1, and is also available on the Github repository (30). The numbers of selected genomes are 2,268, 3,490, 10,587, and 51,676 for the respective collection dates.

MSA.

MSAs were constructed with the online alignment server Multiple sequence alignment using Fast Fourier Transform (MAFFT) (58, 59) for the two smaller datasets with cutoff dates 1 April 2020 and 8 April 2020. To align the two larger datasets with more than 10,000 sequences, a prealigned reference MSA is recommended to accelerate the alignment and reduce the burden on computational resources. Here, we took the collection with cutoff date 8 April 2020 as the prealigned reference MSA for the two largest datasets with cutoff dates 2 May 2020 and 8 August 2020. The MSAs used are given as Dataset S2, and are also available on the Github repository (30).

The MSA is a big matrix S={σin|i=1,…,L,n=1,…,N}, composed of N genomic sequences which are aligned over L positions (16, 21). Each entry σin of matrix S is either one of the four nucleotides (A, C, G, T), or “not known nucleotide” (N), or the alignment gap “-” introduced to treat nucleotide deletions or insertions, or some minorities.

MSA Filtering.

Before filtering, we transform the MSA in two different ways as follows: 1) The gaps “-” are transformed to “N” while the minors “KFY…” are mapped to “N.” Thus five states remain, where “NACG”’ are represented by “12345”; 2) the minors “KFY…” are mapped to “N.” Then there are six states, with “-NACGT” represented by “012345.”

The following criteria are used for prefiltering of the MSA from the 8 August 2020 dataset. If, for one locus, the same nucleotide is found in more than 96.5% of this column, or if the sum of the percentages of A, C, G, and T at this position is less than 20%, then this locus will be deleted. For each sequence, if the percentage of a certain nucleotide is more than 80%, or if the sum of the percentages of A, C, G, and T in this sequence is less than 20%, then this sequence will be deleted. With this filtering criteria, many loci but no sequences are deleted. There are left 51,676 sequences and 689 loci.

B-effective Number.

To mitigate the effects of dependent samplings, it is standard practice to attach to each collected genome sequence σ(b) a weight wb (15, 16, 60), which normalizes its impact on the inference procedure. An efficient way to measure the similarity between two sequences σ(a) and σ(b) is to compute the fraction of identical nucleotides and compare it with a preassigned threshold value x in the range 0≤x≤1. The weight of a sequence σ(b) can be set as wb=1/mb, with mb as the number of sequences in the MSA that are similar to σ(b),mb=|{a∈{1,…,B}}:overlap(σ(a),σ(b))≥x|;[1]here, overlap is the fraction of loci where the two sequences are identical. The B-effective number of the transformed sequences is defined asBeff=∑b=1Bwb.[2]We compare the Beff value with different x for the filtered MSA with q=5 and q=6. As shown in Fig. 5, the dataset with six states shows larger Beff number for all tested x. We thus perform our analysis on the dataset with q=6 states, where “-NACGT” is represented by “012345.”

• Download figure
• Open in new tab
• Download powerpoint

Fig. 5.

Beff number of the 8 August 2020 prefiltered dataset with threshold x. Red denotes q=6 states (“-NACGT”), and black denotes q=5 states (“NACGT”). The number of states is determined by the transform criteria of the prefiltered MSA.

The reweighting procedure partially addresses a point raised (51), that sequenced viral genomes are not a random sample of the global population. That is, even if sequencing is biased by the country they occur in and by contact tracing, sufficiently similar genomes will have lower weight, and so each will contribute less to predictions.

Elements of QLE.

The phenomenon of QLE was discovered by M. Kimura (10) while investigating the steady-state distribution over two biallelic loci evolving under mutation, recombination, and selection, with both additive and epistatic contributions to fitness. In the absence of epistasis, such a system evolves toward linkage equilibrium (LE) where the distribution of alleles at the two loci are independent. The covariance of alleles at the two loci then vanishes. In the presence of pair-wise epistasis and sufficiently high rate of recombination, the steady-state distribution takes form of a Gibbs–Boltzmann formP(σ1,…,σL)=1Zexp{−H(σ1,…,σL)},[3]with an “energy function”H(σ1,…,σL)=∑ihi(σi)+∑ijJij(σi,σj).[4]In the above, Jij can be related to the epistatic contribution to fitness between loci i and j with alleles σi and σj (11⇓–13). The quantity hi is similarly a function of allele σi which depends on both additive and epistatic contributions to fitness involving locus i. It has been verified in in silico testing that, when the terms in Eq. 4 can be recovered, this is a means to infer epistatic fitness from samples of genotypes in a population (14). In the bacterial realm, this approach was used earlier to infer epistatic contributions to fitness in the human pathogens Streptococcus pneumonia (61) and Neisseria gonorrhoeae (62), both of which are characterized by a high rate of recombination. The method was also tested on data on the bacterial pathogen Vibrio parahemolyticus (63). In that study, the results from DCA were not superior to an analysis based on Fisher exact test; see SI Appendix, Different Quantifications of Correlations for a discussion. This is consistent with the approach taken here, as V. parahemolyticus has a low rate of recombination. Further details on the QLE state of evolving populations are given in SI Appendix, Quasi-Linkage Equilibrium (QLE) and Its Range of Validity.

Inference Method for Epistasis between Loci.

The basic assumption of modeling the filtered MSA is that it is composed of independent samples that follow the Gibbs–Boltzmann distribution Eq. 3 with H as in Eq. 4. Higher-order interactions are also possible to include, but we ignore them here (64). This assumption is a simplification of the biological reality; however, it provides an efficient way to extract information from massive data.

On the other hand, in the context of inference from protein sequences, it has been argued that the one encoded in Eqs. 3 and 4 is the minimal generative model, that is, capable not only of reproducing the empirical frequencies and correlations but also of generating new sequences indistinguishable from natural sequences (16, 65, 66).

Many techniques have been developed to infer the direct couplings in Eq. 3, as reviewed in ref. 31 and references therein; see also SI Appendix, Methods of DCA. We employ the PLM method (13, 60, 67⇓⇓–70) to infer the epistatic effects between loci from the aligned MSA. PLM estimates parameters from conditional probabilities of one sequence conditioned on all of the others. For a Potts model with multiple states q>2, this conditional probability isP(σi|σ\i)=exphi(σi)+∑j≠iJij(σi,σj)∑u⁡exphi(u)+∑j≠iJij(u,σj),[5]with u={0,1,2,3,4,5} as the possible state of σi. Eq. 5 depends on a much smaller parameter set compared with that in Eq. 3. This leads to a much faster inference procedure of parameters compared with the maximum likelihood method. With a given independent sample set, one can maximize the corresponding log-likelihood functionPLihi,{Jij}=1N∑shiσi(s)+1N∑s∑j≠iJij(σi(s),σj(s))−1N∑slog∑uexphi(u)+∑j≠iJij(u,σj(s)),[6]where s labels the sequences (samples), from 1 to N. With the filtered MSA, we then run the asymmetric version of PLM (60) in the implementation PLM available in ref. 71 with regularization parameter λ=0.1. The inferred interactions between loci i and j are scored by the Frobenius norm.

Relation to Correlation Analysis.

In LE, the distributions of alleles over different loci are independent. Given unlimited data and unlimited computational resources, the terms Jij in Eq. 4 inferred from the data would then be zero. The locus–locus covariances, defined ascij(a,b)=1σi,a1σj,b−1σi,a1σj,b,[7]would also be zero. The Frobenius norm of cij(a,b) over indices (a,b) as a score of strength of correlations would be zero as well. The qualitative difference between correlation analysis and global model inference based on Eqs. 3 and 4 is that two loci i and j may be correlated (“indirectly coupled”) even if their interaction Jij is zero, provided they both interact with a third locus k. Data in Table 4 and Fig. 4 show that the leading interactions retrieved by DCA cannot be stably recovered in correlation analysis. A different score of statistical dependency between two categorical random variables is mutual information (MI). Results in SI Appendix, Fig. S5 and Table S1 show that the result does not substantially change if using MI instead of Frobenuis norm of correlation matrices. Circos plots of interactions based on correlation scores are available in ref. 30.

Epistasis Analysis with PLM Scores.

PLM procedure yields a fully connected paradigm between pair-wise loci. To extract important information from massive SARS-CoV-2 genomic sequences, we focus on the significant scores between loci, the top 200 pairs. With a reference sequence “Wuhan-Hu-1,” we identify the positions of the corresponding nucleotides. The visualization of these epistases is performed by “circos” software (72).

Randomized Background Distributions.

A way to assess the validity of a small number of leading retained predictions among a much larger set of mostly discarded predictions is to compare to randomized backgrounds. The retained predictions are then, in any case, large (by some measure) and would also be retained if selection were made according to some cutoff, or an empirical p value. The problem is thus how to distinguish the case where a small subset of retained values are large, because they are different, from the case when, in a large number of samples, such values would appear at random. This problem can be addressed by comparing the retained values to the largest values from the same procedure applied to randomized data, as was done for predicted RNA–RNA binding energies in a noncoding RNA discovery pipeline (73). In the context of DCA (PLM) applied to genome-scale MSAs, two earlier studies relying on randomized background distributions are described in refs. 13 and 74.

PLM Scores with Randomization.

To understand the nature of the top 200 PLM scores, we perform two distinct randomization strategies on the MSA, such that its conservation patterns and (or) phylogenetic relations are preserved, while intrinsic coevolutionary couplings (epistatic interactions) are removed (75). Running DCA on artificial sequences ensembles generated by these strategies, and comparing them to the results obtained from original MSA, allows disentangling of spurious couplings given by finite-size effects or by phylogeny. The first strategy, which we refer to as “profile,” randomizes the input MSA by random but independent permutation of all its columns, conserving the single-column statistics for all sites. This destroys all kind of correlations, and DCA couplings inferred from such samples are only nonzero due to the noise caused by finite sample size. In the second strategy, referred to as “phylogeny,” the original MSA is randomized by a simulated annealing schedule where columns and rows are changed simultaneously but so that intersequence distances are kept invariant. Phylogeny inferred from intersequence distance information would therefore be unchanged. Conversely, if the predicted epistatic interactions are due to phylogeny, they should also show up in terms recovered by PLM from MSAs scrambled by “phylogeny.” More details on the randomization strategies can be found in SI Appendix, Phylogenetic Randomization of DCA: Principles.