From sequenced SARS-CoV-2 genomes to a phylogenetic tree
Daniel Lichtblau
Daniel Lichtblau
The Wolfram Data Repository has a curated and daily updated collection of SARS-CoV-2 genome sequences. These sequences are comprised of strings containing the standard nucleotide letters "ACTG". I will show how to create dendrograms from these strings in order to visualize relative genetic proximities both by place and date of collection.
Initial remarks
Initial remarks
This is an updated version of a computational notebook available here. Also it builds on one part of this Wolfram Community notebook from early February.
A related analysis (albeit not with Mathematica as far as I am aware), using similar methods to those I showed, has recently been published:
Gurjit S. Randhawa,, Maximillian P. M. Soltysiak , Hadi El Roz, Camila P. E. de Souza, Kathleen A. Hill, Lila Kari.Machine learning using intrinsic genomic signatures for rapid classification of novel pathogens: COVID-19 case study. PLoSONE 15(4):e0232391, April 24, 2020, 24 pages.
https://doi.org/10.1371/journal.pone.0232391
I will note that I first learned about the Chaos Game Representation, a primary tool in this work, from a 2016 talk given by Lila Kari. Since which time I have used it for species identification, phylogenetic tree construction and, of all things, authorship forensics.
A related analysis (albeit not with Mathematica as far as I am aware), using similar methods to those I showed, has recently been published:
Gurjit S. Randhawa,, Maximillian P. M. Soltysiak , Hadi El Roz, Camila P. E. de Souza, Kathleen A. Hill, Lila Kari.Machine learning using intrinsic genomic signatures for rapid classification of novel pathogens: COVID-19 case study. PLoSONE 15(4):e0232391, April 24, 2020, 24 pages.
https://doi.org/10.1371/journal.pone.0232391
I will note that I first learned about the Chaos Game Representation, a primary tool in this work, from a 2016 talk given by Lila Kari. Since which time I have used it for species identification, phylogenetic tree construction and, of all things, authorship forensics.
Obtaining and processing the sequences
Obtaining and processing the sequences
Our first order of business is to obtain the sequences from the Wolfram Data Repository. They are in the form of a Dataset and can be obtained via this URL. A convenient examples notebook can be found at that link. It has exactly the code I need in order to get lists of sequences along with their collection geographic region and date. Moreover there is a simple way to filter so that we only use fully sequenced examples. This code, and the data resource, are the work of my colleague John Cassel.
Obtain the sequences and convert to a Mathematica List structure. This download is from April 30, 2020.
fullGenes=ResourceData["Genetic Sequences for the SARS-CoV-2 Coronavirus"][Select[StringContainsQ[#GenBankTitle,"complete genome"]&],{"Sequence","GeographicLocation","CollectionDate"}];geneLists=Normal[fullGenes];Length[geneLists]
Out[]=
1389
As it happens, I already know that many have a bunch of trailing adenine nucleotides, and the presence or absence thereof does not correspond to mutations. The resource notebook has convenient code for removing these:
In[]:=
dropTrailingA[seq_]:=StringReplace[seq,StartOfString~~Shortest[a__]~~("A"..)~~EndOfStringa];
It is useful to check the size range of the sequence strings:
ReverseSortBy[Tally[Map[StringLength[dropTrailingA[#]]&,"Sequence"/.geneLists]],First]
Out[]=
{{29923,1},{29917,1},{29914,1},{29912,2},{29910,1},{29909,2},{29907,1},{29906,1},{29904,1},{29903,1},{29902,1},{29901,2},{29899,1},{29897,2},{29894,2},{29891,1},{29890,1},{29888,1},{29887,1},{29883,2},{29882,1},{29881,5},{29880,1},{29879,3},{29877,3},{29876,3},{29875,1},{29873,2},{29872,2},{29871,2},{29870,278},{29869,15},{29868,25},{29867,21},{29866,18},{29865,7},{29864,31},{29863,17},{29862,11},{29861,9},{29860,20},{29859,20},{29858,14},{29857,17},{29856,15},{29855,7},{29854,12},{29853,9},{29852,11},{29851,10},{29850,6},{29849,5},{29848,8},{29847,7},{29846,16},{29845,8},{29844,8},{29843,9},{29842,11},{29841,6},{29840,15},{29839,66},{29838,18},{29837,6},{29836,14},{29835,22},{29834,17},{29833,21},{29832,13},{29831,21},{29830,13},{29829,9},{29828,3},{29827,3},{29826,3},{29825,5},{29823,6},{29822,2},{29821,9},{29820,4},{29819,7},{29818,2},{29817,4},{29816,2},{29815,5},{29814,4},{29813,9},{29811,4},{29810,5},{29808,4},{29807,2},{29806,1},{29805,1},{29803,2},{29802,1},{29801,1},{29800,1},{29799,2},{29798,2},{29797,5},{29796,10},{29795,1},{29793,3},{29792,1},{29791,1},{29790,7},{29789,2},{29788,5},{29787,3},{29786,1},{29785,4},{29784,1},{29783,1},{29782,16},{29781,15},{29780,2},{29779,1},{29778,1},{29777,1},{29776,6},{29775,1},{29774,3},{29772,1},{29771,4},{29768,1},{29767,1},{29764,2},{29763,1},{29761,3},{29760,2},{29759,2},{29757,1},{29756,12},{29755,1},{29754,1},{29750,1},{29749,1},{29748,1},{29747,2},{29746,3},{29745,2},{29744,6},{29743,1},{29740,2},{29739,2},{29738,2},{29737,1},{29736,4},{29732,4},{29731,3},{29730,3},{29729,6},{29727,3},{29726,22},{29725,2},{29724,44},{29723,9},{29722,5},{29720,4},{29719,5},{29718,2},{29717,4},{29716,6},{29715,6},{29714,8},{29713,3},{29712,4},{29711,1},{29710,2},{29709,1},{29707,1},{29706,2},{29705,4},{29704,3},{29702,1},{29701,1},{29700,2},{29697,1},{29696,1},{29693,1},{29688,2},{29687,2},{29685,1},{29683,1},{29681,1},{29679,3},{29678,1},{29673,1},{29671,2},{29664,1},{29662,2},{29658,1},{29651,1},{29649,1},{29643,1},{29629,1},{29619,1},{29618,1},{29617,1},{29615,1},{29613,3},{29611,5},{29610,2},{29608,2},{29596,1},{29589,1},{29583,1},{29574,1},{29554,1},{29488,1},{29432,1}}
Removing outliers
Removing outliers
Of these many genome sequences, 278 have the same length of 29870 (adjusted for dropped trailing A sequences). It is easy to see that most are in a fairly tight range. We will initially restrict to these similar length sequences on grounds that the sequencing is not known to be perfect and we expect the lesser-quality sequenced examples to fall outside the size range. Also I was informed by my colleague John Cassel that the lone Netherlands sequence was taken from a European mink, not a human. So we remove that one as well.
Obtain lists with entries of the form (sequence, location, date):
In[]:=
normalGL=Map[Values,geneLists];
Prune off those with sizes outside the range where most are located and see how many remain.
midsize0=Select[normalGL,29830<=StringLength[dropTrailingA[#[[1]]]]≤29880&];midsize=Select[midsize0,FreeQ[#,"Netherlands"]&];Length[midsize]
Out[]=
893
Obtain the sequences, locations, and dates as three separate lists:
In[]:=
{fullseqs0,places,dates}=Transpose[midsize];fullseqs=Map[dropTrailingA,fullseqs0];
Tally the sequence lengths for the retained genomes:
ReverseSortBy[Tally[Map[StringLength,fullseqs]],First]
Out[]=
{{29880,1},{29879,3},{29877,3},{29876,3},{29875,1},{29873,2},{29872,2},{29871,2},{29870,278},{29869,15},{29868,25},{29867,21},{29866,18},{29865,7},{29864,31},{29863,17},{29862,11},{29861,9},{29860,19},{29859,20},{29858,14},{29857,17},{29856,15},{29855,7},{29854,12},{29853,9},{29852,11},{29851,10},{29850,6},{29849,5},{29848,8},{29847,7},{29846,16},{29845,8},{29844,8},{29843,9},{29842,11},{29841,6},{29840,15},{29839,66},{29838,18},{29837,6},{29836,14},{29835,22},{29834,17},{29833,21},{29832,13},{29831,21},{29830,13}}
Creating a color list to correspond to collection data
Creating a color list to correspond to collection data
Grouping by collection date
Grouping by collection date
We will create ordinary numbers from the collection dates. Since this is only a crude indicator of when a virus was acquired we do nothing fancy here. The collection times are given to the hour but we only use the day. Since these first appear in 2019 we subtract that from the year and multiply by 365, then multiply the month number by the average number of days (I use 30.4 for this), then add the day of month.
First see how a date is formed:
dates[[1]]//InputForm
Out[]//InputForm=
DateObject[{2020, 3, 6}, Day, Gregorian, -5.]
So we take the first element from each of these, perform the above-mentioned arithmetic, and then rescale so that the numbers range from 0 to 1:
In[]:=
dates1=Map[First,dates];ndates1=Map[(#[[1]]-2019)*365+#[[2]]*30.4+If[Length[#]3,#[[3]],15]&,dates1];ndates=Rescale[ndates1]+.02;
The rescaling was so that we could conveniently assign colors to dates:
colorDateLabels=Thread[Style[dates1,Map[Darker[Hue[#]]&,ndates]]];Short[colorDateLabels]
Out[]//Short=
{{2020,3,6},{2020,3,31},{2020,3,31},{2020,3,24},{2020,3,24},{2020,3,31},882,{2020,1,19},{2020,1,11},{2020,1,10},{2019,12},{2019,12}}
Grouping by collection location
Grouping by collection location
From inspection I know that every collection location is either a geographic Entity or else of the form Missing["Unavailable"]. So it is easy to extract the location names to use as dendrogram node labels. Of course many infection samples collected in one location may have originated elsewhere. And some geographic locations are not well "localized" e.g. China or United States. So we hardly expect perfect correlation between genome proximity and collection location. But we do expect that there will be some visible evidence of a relation between genetic and geographic similarity.
placeNames=Map[If[Head[#]===Entity,#[[2]],Head[#]]&,places];placeSet=Union[placeNames]
Out[]=
{Australia,Brazil,China,France,Greece,HongKong,India,Israel,Italy,Japan,Malaysia,Pakistan,Peru,SouthAfrica,SouthKorea,Spain,SriLanka,Sweden,Taiwan,Turkey,UnitedStates,Vietnam}
Again we form a color scheme:
In[]:=
nplaces=Length[placeSet];crange=Rescale[Range[nplaces]/(nplaces+1)]+.02;placeColors=Map[Darker[Hue[#]]&,crange];
Now associate colors to geographic locations:
colorPlaceRule=Thread[placeSet->placeColors];colorPlaceLabels=Thread[Style[placeNames,placeNames/.colorPlaceRule]];Short[colorPlaceLabels]
Out[]//Short=
{China,UnitedStates,UnitedStates,UnitedStates,UnitedStates,UnitedStates,882,UnitedStates,China,China,China,China}
Put these together.
colorPlaceDateLabels=Thread[{colorPlaceLabels,colorDateLabels}];Short[colorPlaceDateLabels]
Out[]//Short=
{{China,{2020,3,6}},{UnitedStates,{2020,3,31}},{UnitedStates,{2020,3,31}},888,{China,{2019,12}},{China,{2019,12}}}
Creating a dendrogram with labels given by collection date
Creating a dendrogram with labels given by collection date
A phylogenetic tree is a type of dendrogram one can use to visualize proximities and branching of gene sequences. In order to create such a graphic there must be an underlying way to obtain "distances" between pairs of sequences. One approach would be to use string distance measures. This can be quite slow when the strings have close to 30,000 letters. I developed a surrogate based on a well-known method for turning gene sequences into numeric arrays along with some standard methods for dimension reduction. I will defer on the details until a later section since we will actually need to make alterations. For now I simply use the phylogenetic tree plotting Wolfram Resource Function "out of the box". The tree does seem to show clumps grouped by collection date.
AbsoluteTiming[dgram=ResourceFunction["PhylogeneticTreePlot"][fullseqs,colorPlaceDateLabels];]
Out[]=
{147.659,Null}
Show[dgram,AspectRatio14,ImageSize800]
Out[]=
A few things should be noted. One is that some of the positioning is heavily tied to the distance function used in the clustering. Another is that there are a number of different methods that can be used in clustering. One can use the Method and/or DistanceFunction options of PhylogeneticTreePlot to obtain different results. Also some subtrees might appear to be close due to contiguity in the plot, when a rotation would make them appear far apart, or vice versa. Finally there is a disclaimer: We use geographic administrative divisions as provided by the data sources, and make no claim as to nationalities or borders. I follow that convention in the remarks below.
While there may be a few surprises, grouping is overall not unexpected. I will mention some groupings I find interesting, starting from the top.
1
.Two sequences from China, in December and January, are in the midst of US genomes from March and April.
2
.There is a small cluster from Turkey, Sri Lanka, Taiwan and Hong Kong. Those from Hong Kong clearly precede the others by date.
3
.Sequences from India, France, Sri Lanka, Taiwan, Spain and South Africa, all from March/April, are amid slightly earlier sequences from the US. It would be interesting to know if they got infected in the States.
4
.A pair of sequences from Peru and Israel, both in March, are close to one another and to some US sequences from April. Perhaps notable is that a sizable group of Israelis were temporarily stranded in Peru for some days when flights were stopped.
5
.There is then a fairly long subtree almost entirely from the US, with a handful interspersed from Taiwan. There seem to be clusters based on similar dates. It would be interesting to see how these break down in terms of individual states. This is something I need to check.
6
.There appears to be a subtree mostly from March comprised of several genomes from Greece and Taiwan, as well as some from Spain and Sri Lanka, all amidst many from the States. Several of these last were on the earlier side so it would again be interesting to know if travel from the US was involved.
7
.A subtree from the US contains a few sequences from Hong Kong, China and India. They precede most of the ones from the US by date.
8
. A small subtree consists mostly of US sequences, primarily from March. It contains a handful of late January sequences from Hong Kong, China, Japan, Greece and Taiwan. With the exception of the one from Greece in mid-March, these were all early sequences, dating from January.
9
.A cluster with several from the US and elsewhere in March also contains older sequences from South Korea, Australia, Italy, Sweden and Hong Kong.
10
.A bottom-most subtree is dominated by early sequences, from December, January, and early February. Many of these are from China along with Taiwan, India, Vietnam, the US, and Hong Kong, with a number of others appearing in the US and a few other places from mid-February to mid-March (there are also a few later ones in this subtree).
Description of the underlying method
Description of the underlying method
A happenstance failure in the earlier Computational Essay pointed me in a direction that necessitated a modification to the code in that Wolfram FunctionRepository entry. I will give a brief tour of the inner workings and then explain this modification.
Given a genetic sequence, a method to convert to an image was presented by Joel Jeffrey in the late 80's. A modest refinement, which I use, is called the Frequency Chaos Game Representation (FCGR). We have a convenient Wolfram Function Repository item for this, FCGRImage. I will illustrate it on the first and last of the sequences in our set of SARS-CoV-2 genomes.
A next step is to reduce dimension using the Discrete Cosine Transform (DCT). The images above have dimensions 128x128 and I reduce to 30x30. The idea is to retain only the lower frequency components, in effect blurring the images by removing finer details.
Compute the DCTs:
In order to visualize the effect of this reduction we must invert these (dimension-reduced) transforms:
One will notice the coarse similarities between original FGCR images and those that are dimension-reduced by DCT. The point though is not to visualize these but rather to utilize them in a further dimension reduction. We now process all the genome sequences as above and flatten each into a vector. This gives an array of dimensions 96x900 (since 30x30 is 900).
We now use the Singular Value Decomposition (SVD) to reduce, this time to dimension 40. Here is a helper function I grabbed and slightly modified from PhylogeneticTreePlot implementation code:
Compute the next step in dimension reduction:
Recall that we began with genome strings containing nearly 30K nucleotides (for which, as we all know, the correct notation is 30000 bpi). We then made them into 128x128 pixel images. If you work it out properly, this is far from is a size reduction. We next reduced to 30x30 and then to 40 machine doubles (each taking 8 bytes). This is a substantial reduction. The point of the SVD is to concentrate most of the "energy" (properly defined) of the input into a smaller dimension. With a bit of work we can visualize these reduced vectors.
Here are the pictures; don't sweat the code details:
The important take-away is that these are visually quite similar to those DCT images. Actually there is an important difference I did not yet mention. What I stumbled across earlier is that the SVD packs so much oomph (apologies for the technical term) into the first component that one swamps the rest. To see this, I retain the singular values; here are the initial (largest) of the singular values from the SVD.
That first one is around 75 times larger than the next. What this means for our purposes is this. If we retain it (as is done in PhylogeneticTreePlot) then those dimension-40 vectors would all have large first components and piddly (another technical term, again sorry) remaining values. Since the dendrogram code works by clustering based on distances, we really do not get a reasonable sense of separation because the first components are all about the same and all swamp the rest. So this, finally, explains why I felt a need to modify the inner workings of that resource function. So I added an option to PhylogeneticTreePlot, with an Automatic default setting, to allow removing first components when the first singular value is much larger than the next. I should mention that this issue of a very large first singular value does not tend to arise when working with a set of genomes from different species. Which was the case for most of the experiments I did in a "spare time" project from a year ago and hence why I did not consider it before.
A dimension-reduced point plot
A dimension-reduced point plot
We show the point plot on the vectors formed without that large singular value, using MultidimensionalScaling from the Wolfram Function Repository. Colors are blended from collection dates and locations.
Here are some different views of the points.
One will observe that this plot does not show very distinct clusters, though it might indicate four or more “lineages”. And it does tend to show that genomes tend to have neighbors with similar collection dates and locations.
Summary
Summary
Three-plus months ago I began playing with genome sequences that were being made available for SARS-CoV-2 (which at that time was then still informally known as 2019-nCoV). I took the few sequences then available in the GenBank repository at NCBI, also signed up on GISAID and found several there, and moreover used many other coronavirus genomes from GenBank, for purposes of performing comparisons and various other computations. Much of what resulted can be found in a Wolfram Community post from early February. At that time, one interest of mine was to show how a phylogenetic tree could be used to indicate which coronaviruses were most closely related to the novel SARS-CoV-2. While such analysis had already been done elsewhere, it was interesting to see that methods available in the Wolfram language were suitable for independently (and speedily) obtaining similar results. And such results can give useful clues to possible origins of this virus. For example, I recall reading a post in late January that mentioned "bats, not snakes" (a phrase that has since been repeated elsewhere). A phylogenetic tree containing coronavirus genomes from bats, snakes, and several other species showed this quite clearly. In this note I tackle a problem related to, but more difficult than, the phylogenetic tree analysis in the prior work. The goal herein is to more accurately compute dendrograms for genome samples from the same virus, specifically the SARS-CoV-2 full sequences now available in the Wolfram Data Repository.
This has been an interesting investigation. I am not sure I learned anything new about SARS-CoV-2, but as a lay person I did not really expect to. I did learn something about PhylogeneticTreePlot (an area more familiar to me) that was not really obvious in experiments on genome sets with greater variability than the one we looked at here. This lead us into a brief tour through a method for converting genome sequences into vectors. It has applicability to genome look-up and proximity measures, as well as other aspects of alignment-free genomic computation. And it is an approach to which I am quite partial. I should be, since I spent significant time to develop it over the past year or two. While I'd rather we did not have a test set of quite this significance, it seems to have use for gauging distances in SARS-CoV-2 sequences.
Related Links
Related Links
◼
A walk-through of the SARS-CoV-2 nucleotide Wolfram resource by John Cassel
A parting dendrogram
A parting dendrogram
For those who might be interested, I include the dendrogram created from the full set of fully sequenced genomes. It requires the better part of five minutes of computational time. That could be reduced by using parallel kernels to create the FCGR images and extract the low dimensional frequences from their DCTs.