SARS-CoV-2 in Australia during 2020
SARS-CoV-2 in Australia during 2020
Early in this year one of us (DL) wrote up work on use of alignment-free methods in genomics to study various aspects of the SARS-CoV-2 virus. One example was to plot a three-dimensional representation that indicates changes in the genome over time. One region chosen for this purpose was Australia. The genomes had been collected between February and October of 2020, with the bulk having been from July and August. Here is the picture that resulted.
One will notice a spread over time, in multiple directions. Such trends might indicate genetic mutations that are relatively benign, though of course we do know that certain mutations, in various parts of the world, have given rise to more genetically fit variants. We decided to revisit Australia to check what are the trends that have appeared in 2021. This was motivated by several considerations. One is that the country is geographically isolated in terms of border controls. Another is that the population is not very large, and moreover the percentage of genomes sequenced from patients with COVID-19 is relatively high, so we get a solid data set without being swamped by the data (our method can handle a few thousand data points, but not tens of thousands without making substantial adjustments). A third is that the country does have some contact with the rest of the world, and the variant of Concern have appeared there. Delta, in particular, is now fairly widespread in at least one state, New South Wales (where Sydney is located).
The idea of this study was to get some sense of how widespread these variants are in Australia. What we found, however, has interest beyond that: there may be a new COVID-19 subvariant emerging. Before proceeding to the details we remark that the initial analysis, using alignment-free methods, was also done using recent sequences from Ireland and from Finland. Those did not show signs of a separate variant. It is important also to state that new variants do not as a rule equate to variants of concern. New ones no doubt emerge simply due to random mutations that might spread briefly (think super-spreader event). Unless the mutation set confers advantage to the virus, these can just as quickly die out as more virulent strains (think Delta) take over.
Brief remark on nomenclature: There is no universal agreement amongst news media as to use of capitalization for the variants. That applies to the authors of this post as well.
The four Variants of Concern, plus one relative
The four Variants of Concern, plus one relative
Before showing how Australian sequences fit in, we first indicate using both Multidimensional scaling (MDS) and a phylogenetic tree plot that the methods used are sufficiently discriminating to actually separate the variants from one another. This is important since most alignment-free methods used in the COVID-19 literature to date do not go this far but rather only try to separate SARS-CoV-2 from other coronavirus types. As the methodology behind these two plots has been covered in two articles (this one as well as the one linked to above), a Wolfram blog, and several posts on the Wolfram Community forum, we will not repeat it here.
For our purposes we use recently collected genomes from around the world that have been placed in the GISAID repository. We selected sequences from the latter part of July, going back as far as necessary in order to obtain between 400 and 500 sequences for each of the variants Alpha, Beta and Gamma. Due to the widespread reporting of Delta genomes we had to use over 600 of those; a day later would have resulted in too few for our needs. We mention that we also created similar plots using several vVariants of Interest (specifically, Epsilon, Iota and Lambda). Results indicated these have no presence in Australia so we omitted them from further consideration. But...we learned that VoI Kappa is a close relative of Delta and is present in Australia. So we include it as well.
Here is the 3-D MDS picture. The variants do appear to be well separated even after dimension reduction, except for some intermingling at the edge between sibling variants Delta and Kappa.
We decimated by a factor of four and squeezed the resulting tree. It is again clear that the variants are all nicely separated. Though strangely enough, Kappa seems to be closer to Alpha than it is to Delta ( which might simply indicate a deficiency in how the relative distances are gauged).
Australia 2021
Australia 2021
We now add the GISAID sequences from Australia that were collected starting January 1 and submitted by July 30. There were over 1200 such. We color code by date for these, with specific colors as above, that do not overlap those used for the variants.
There are several things to notice here. One is that many Australian sequences from roughly March and April seem to cluster among the Alpha and Beta variants. Many other older sequences are apart from the clusters. The Delta grouping, which has the largest spread in this picture, seems to cover most of the newer Australian sequences. But there is a recent grouping near the top and to the right of center that seems to be apart from Delta. Another recent clump, near bottom center, lies slightly off the pink Kappa region. (Full and slightly embarrassing disclosure: originally the second author did this without using Kappa samples, and the cluster near bottom center was then standing alone. This is what gave us reason to take the analysis further. Only later, when we learned the Kappa strain was present in Melbourne, did we recognize the need to include it in this study.)
As with all dimension reduction schemes, MDS can be misleading insofar as it in effect projects onto a lower dimension. We need to check a phylogenetic tree in order to get a better idea of the extent to which these clusters might stand apart. We now use color to indicate date, and mark each Australian genome by the state in which it was collected (some of these have further information which we can ignore here). Now we decimate each group by a factor of four, just to keep size and speed within reason.
There is a distinct branch at the top, containing, in addition to the Kappa samples, many recent sequences primarily collected in the state of Victoria (home to the city of Melbourne). Now Victoria has had a modest uptick in reported COVID-19 cases going back to June, and some are claimed to be from Kappa. The alignment-free methods we have used indicate a distinct clustering of the Victoria subset that coincides heavily with both time and geography, but this by itself can only suggest the possibility of a (sub-)variant emerging. We need finer analysis to see more specifically what is happening in the genomes within this cluster. For this purpose we require slower but more fine-tuned alignment-based methods.
Aligning SARS-CoV-2 sequences
Aligning SARS-CoV-2 sequences
Sequence alignments provide the exact differences between genetic sequences. By aligning a set of sequences against a reference sequence, we can then see exactly what differences cause them to be more or less similar to each other. We knew that in order to understand this new cluster in Victoria, we had to see what genetic differences were in common between members of the cluster, and if those were different than other variants.
Despite mutation, the genetic sequences of any complete sample of SARS-CoV-2 and the SARS-CoV-2 reference sequence are nearly the same. Fortunately, this means it is quite practical to gather the differences between a sample and the reference, as various optimizations can be applied to find how the sequences fit together, and each sample can ultimately be represented by a their variations alone. As it happens, we had the tools ready in the Wolfram Function Repository to look into these differences, as will be describe soon in a forthcoming blog. Said briefly, for all of the sample sequences under study, we aligned them against the reference sequence, extracted differences from those alignments, and removed those differences explained by lack of certainty in the sample data.
For each sample, this process produces a list of differences, where each difference is of the form “position:reference>variation”. As a made-up example, “534:A>”C” would say that at position 534 in the reference sequence, substitute “C” for “A” to obtain the variant sequence. Insertions and deletions can be expressed as empty strings on either the left or right side of the rule, respectively. This is a quite compressed representation, where an example delta variant sample of length 29851 can be represented with only 42 differences, leading to it only requiring 6% of the space compared to the complete string.
For three sequences in the cluster, chosen to most closely embody the cluster, we found the differences they all had in common, and then removed differences found in example sequences from the Alpha, Delta, and Kappa variants. This not only allows us to rule out the cluster being a result of those variants, but also filter out commonly-occurring background variations that developed after the reference sequence was formalized but are now widespread.
We then looked to the number of the variations that were found in each of the Australian sequences, and our results we very interesting. When more than a handful of variations were found, they were nearly all found.
We then looked at labels corresponding to the the difference lists identified by these variations, and found they were indeed found in the cluster. We then gathered all of the variations in the difference lists found in the cluster identified by this method, ordered them by total proportion found, and took the most frequent to establish the overall set of cluster variations.
Given that we had found a set of variations that clearly identified the cluster, we were interested to see how they compared against other variations. We found that it shared no frequent variations with alpha, beta, or gamma. However, the cluster did share a fair number of variations with both delta and kappa, though the cluster had some variations found in neither. It seems that most to all of the variations found both in this cluster and delta are also found in kappa.
Cluster mutations not found in Kappa are:
Kappa mutations not found in the cluster:
We also looked into whether or not the variations of the unknown cluster were in protein coding regions of the viral genome, and if they would make a coding difference. Indeed, a reasonably substantial number of them were coding variations, including quite a few on the surface (“spike”) protein, which is arguably the most relevant to immune response. When we removed the coding variations found in delta and kappa, there were only a few coding variations, although one of them was a surface variation.
The amino acid substitutions specific to the cluster from Victoria are given in the table below.
Conversely, substitutions common to Kappa but not the Victoria cluster are in the next table.
It should be noted that one change in the cluster, and several in Kappa, involve the surface glycoprotein, better known as the spike protein. Such difference can be important in terms of how they affect fitness of the virus e.g. via transmissibility or pathogenicity. Also the change specific to the Australia cluster, H1101Y, is common to the P.3 lineage from the Philippines. See for example this article or this other (where the variant is referred to as PHL-B.1.1.28) for more details (the confusion in nomenclature is explained by the Wikipedia article on the Theta variant). In contrast, the Kappa variant instead shows the substitution H1101D in about 22.3% of the sequences used in this study. The I33T mutation has also been associated with severe forms of the illness. It is apparently common to the B .1.1.33 lineage.
Concluding remarks
Concluding remarks
This brief (except for pictures!) essay shows ways in which we can utilize Mathematica to analyze genome sequences. These go from broad-based comparison all the way down to the level of specific nucleotide mutations and corresponding amino acid substitutions. We apply these to genomes responsible for the current COVID-19 pandemic. What we found was something of a surprise, insofar as a cluster of cases from Victoria, Australia, which were claimed to be of the Kappa variant, in fact turn out to have mutations distinct from Kappa and to lack some of the mutations common to Kappa. After some polishing we hope to make the full notebook available in the near future (albeit sans the actual genomes as these are available only from GISAID).
On a related note, an upcoming blog (same authors) will illustrate use of material in the Wolfram Function Repository for purposes of genomic analysis. It will cover much of the methods used in this post.
GISAID acknowledgements
GISAID acknowledgements
Elbe, S., and Buckland-Merrett, G. (2017) Data, disease and diplomacy: GISAID’s innovative contribution to global health. Global Challenges, 1:33-46. DOI:10.1002/gch2.1018 PMCID: 31565258
Shu, Y., McCauley, J. (2017) GISAID: Global initiative on sharing all influenza data – from vision to reality. EuroSurveillance, 22(13) DOI:10.2807/1560-7917.ES.2017.22.13.30494 PMCID: PMC5388101