Collecting data

Accession list of SARS-like viruses:
In[]:=
viruses={"AY274119.3","AY278488.2","AY278491.2","AY278554.2","AY278741.1","AY282752.2","AY283794.1","AY283795.1","AY283796.1","AY283797.1","AY283798.2","AY291451.1","AY304488.1","AY390556.1","AY394985.1","AY394996.1","AY395003.1","AY502927.1","AY559093.1","EU371559.1","EU371562.1","EU371564.1","GU553365.1","KF294455.1","KJ473816.1","KY417144.1","KY417145.1","KY967357.1","KY983587.1","MF542265.1","MG772933.1","MG772934.1","MK211376.1","NC_002645.1","NC_004718.3","NC_005831.2","NC_019843.3","NC_038294.1","MN996532","KF367457"};
One can use the function ImportFASTA, from Wolfram Function Repository, to Import FASTA data from the NCBI using an NCBI Reference Sequence. To save time, I have done it and gave it an iconized form.
In[]:=
virusNameSeq=
VirusNameSequece
(*Map[ResourceFunction["ImportFASTA"],viruses]*);
Take 1st element of from all sublist of virusNameSeq, flatten the result, and return the substrings from the first space to “,”:
In[]:=
names=Flatten[StringCases[virusNameSeq[[All,1]]//Flatten," "~~x__~~","x]]
Out[]=
{Severe acute respiratory syndrome-related coronavirus isolate Tor2,SARS coronavirus BJ01,SARS coronavirus HKU-39849,SARS coronavirus CUHK-W1,SARS coronavirus Urbani,SARS coronavirus CUHK-Su10,SARS coronavirus Sin2500,SARS coronavirus Sin2677,SARS coronavirus Sin2679,SARS coronavirus Sin2748,SARS coronavirus Sin2774,SARS coronavirus TW1,SARS coronavirus SZ16,SARS coronavirus GZ02,SARS coronavirus HSZ-Bb,SARS coronavirus ZS-B,SARS coronavirus ZS-C,SARS coronavirus TW4,SARS coronavirus Sin845,SARS coronavirus ZJ02,SARS coronavirus BJ182-4,SARS coronavirus BJ182-12,SARS coronavirus HKU-39849 isolate TCVSP-HARROD-00003,SARS-related bat coronavirus isolate Anlong-111 orf1ab polyprotein and orf1a polyprotein genes,BtRs-BetaCoV/YN2013,Bat SARS-like coronavirus isolate Rs4084,Bat SARS-like coronavirus isolate Rf4092,Human coronavirus 229E strain SC2872,Human coronavirus 229E strain SC3112,Human coronavirus 229E strain 229E/Haiti-1/2016,Bat SARS-like coronavirus isolate bat-SL-CoVZC45,Bat SARS-like coronavirus isolate bat-SL-CoVZXC21,Coronavirus BtRs-BetaCoV/YN2018B,Human coronavirus 229E,SARS coronavirus,Human Coronavirus NL63,Middle East respiratory syndrome coronavirus,Betacoronavirus England 1,Bat coronavirus isolate RaTG13,Bat SARS-like coronavirus WIV1}
Replacement rule to change U by T and keep only A, C, G and T:
In[]:=
srules={"U""T",Except[Characters["ACGT"]]""};
Give an association with viruses accession as keys and sequences as values, then replace the sequences using srules:
In[]:=
dataSARSCoV=StringReplace[srules]/@AssociationThread[viruses,Flatten[virusNameSeq[[All,2]]]];
Take COVID-19 data from WDR, then select those whole GenBankTitle contains “complete genome””
In[]:=
dataWDR=ResourceData["Genetic Sequences for Novel Coronavirus 2019-nCoV from Wuhan, China"][Select[StringContainsQ[#GenBankTitle,"complete genome"]&]];
From previous imported data, take only Accession and Sequence columns, make it as association with accession as keys and sequence as values, finally replace the sequences using srules:
In[]:=
dataCOVID19=StringReplace[srules]/@<|Rule@@@Normal[Values[dataWDR[All,{"Accession","Sequence"}]]]|>;
Generate an association from SARS-like and COVID-19 viruses with accessions as keys and corresponding names as values:
In[]:=
accessionNameAssoc=Join[AssociationThread[viruses(Style[#,Bold,Italic]&/@names)],AssociationThread[Keys[dataCOVID19]StringSplit[Normal[dataWDR[All,"GenBankTitle"]],","][[All,1]]]];

Z-curve

For a given genetic sequence of length N, let’s denote the number of occurrence of A, C, G, and T bases at the n-th subsequence by
A
n
,
C
n
,
G
n
, and
T
n
. Then, the corresponding 3D coordinate for the n-th subsequence is obtained by the following transformation:
Out[]//TraditionalForm=
x
n
y
n
z
n

1
-1
1
-1
1
1
-1
-1
1
-1
-1
1
A
n
C
n
G
n
T
n
For example, for a sequence such as AACGT, we have:
{
A
1
=1,
C
1
=0,
G
1
=0,
T
1
=0}
,
{
A
2
=2,
C
2
=0,
G
2
=0,
T
2
=0}
,
{
A
3
=2,
C
3
=1,
G
3
=0,
T
3
=0}
,
{
A
4
=2,
C
4
=1,
G
4
=1,
T
1
=0}
, and
{
A
5
=2,
C
5
=1,
G
5
=1,
T
5
=1}
, and the corresponding Z-curve would be {{0,0,0},{1,1,1},{2,2,2},{1,3,1},{2,2,0},{1,1,1}}.
Note that since
x
n
=(
A
n
+
G
n
)-(
C
n
+
T
n
)
it is the different between the Purine bases (R) and Pyrimidine ones (Y). Therefore, it is a measure of RY disparity. Likewise,
y
n
quantifies the Amino (M) and Keto (K) difference, called as MK-disparity and
z
n
the WS-disparity (W for Weak and S for Strong Hydrogen bonds).
Generate Z-curve by cumulative sums of four matrices assigned to A, C, G and T:
In[]:=
ZcurveGenerator[st_]:={{1,-1,1,-1},{1,1,-1,-1},{1,-1,-1,1}}.#&/@FoldList[#1+Switch[#2,"A",{1,0,0,0},"C",{0,1,0,0},"G",{0,0,1,0},"T",{0,0,0,1}]&,{0,0,0,0},Characters[st]]
Generate Z-curves of COVID-19 viruses (the result is an association in the form of accession->Z-curve):
In[]:=
zCurveCOVID19=ZcurveGenerator/@dataCOVID19;
Generate Z-curves of SARS-like viruses (the result is an association in the form of accession->Z-curve):
In[]:=
zCurveSARS=ZcurveGenerator/@dataSARSCoV;
Represent 3D image of Z-curve where the red ones belongs to SARS-like viruses and blue ones of COVID-19:
Legended[Graphics3D[{Blue,Values[Line/@zCurveCOVID19],Red,Values[Line/@zCurveSARS]},AxesTrue,ImageSizeLarge,BoxRatios{2,4,4},AxesLabel{"RY disparity","MK disparity","WS disparity"}],LineLegend[{Blue,Red},{"Z-curves of 39 COVID19 viruses","Z-curves of 40 SARS-like viruses"}]]
Out[]=
Find the shortest sequence and denote its length by “min”:
In[]:=
min=Min[Min[Length/@zCurveCOVID19],Min[Length/@zCurveSARS]];
Trim the Z-curves up to the length of shortest sequence:
In[]:=
trimmedZCurveCOVID19=Take[#,min]&/@zCurveCOVID19;​​trimmedzCurveSARS=Take[#,min]&/@zCurveSARS;

Clustering and classification based on Z-curves

Find the 3 nearest Z-curves of SARS-like viruses to the Z-curve of each COVID-19 viruses:
In[]:=
closest=Nearest[trimmedzCurveSARS,Values[trimmedZCurveCOVID19],3];
In[]:=
Grid[Prepend[Flatten/@Transpose[{Keys[trimmedZCurveCOVID19],closest}],{Style["COVID-19 Accession",Bold],Style["3 nearest SARS-like viruses",Bold],SpanFromLeft,SpanFromLeft}],FrameAll]
Out[]=
COVID-19 Accession
3 nearest SARS-like viruses
MN908947
MN996532
MG772934.1
MG772933.1
NC_045512
MN996532
MG772934.1
MG772933.1
MN938384
MN996532
MG772934.1
MG772933.1
MN975262
MN996532
MG772934.1
MG772933.1
MN985325
MN996532
MG772934.1
MG772933.1
MN988713
MN996532
MG772934.1
MG772933.1
MN988668
MN996532
MG772934.1
MG772933.1
MN988669
MN996532
MG772934.1
MG772933.1
MN994467
MN996532
MG772934.1
MG772933.1
MN994468
MN996532
MG772934.1
MG772933.1
MN997409
MN996532
MG772934.1
MG772933.1
MN996527
MN996532
MG772934.1
MG772933.1
MN996528
MN996532
MG772934.1
MG772933.1
MN996529
MN996532
MG772934.1
MG772933.1
MN996530
MN996532
MG772934.1
MG772933.1
MN996531
MN996532
MG772934.1
MG772933.1
MT007544
MN996532
MG772934.1
MG772933.1
MT019529
MN996532
MG772934.1
MG772933.1
MT019530
MN996532
MG772934.1
MG772933.1
MT019531
MN996532
MG772934.1
MG772933.1
MT019532
MN996532
MG772934.1
MG772933.1
MT019533
MN996532
MG772934.1
MG772933.1
MT020880
MN996532
MG772934.1
MG772933.1
MT020881
MN996532
MG772934.1
MG772933.1
LC521925
MN996532
MG772934.1
MG772933.1
MT027062
MN996532
MG772934.1
MG772933.1
MT027063
MN996532
MG772934.1
MG772933.1
MT027064
MN996532
MG772934.1
MG772933.1
LC522972
MN996532
MG772934.1
MG772933.1
LC522973
MN996532
MG772934.1
MG772933.1
LC522974
MN996532
MG772934.1
MG772933.1
LC522975
MN996532
MG772934.1
MG772933.1
MT039873
MN996532
MG772934.1
MG772933.1
MT039887
MN996532
MG772934.1
MG772933.1
MT039888
MN996532
MG772934.1
MG772933.1
MT039890
MN996532
MG772934.1
MG772933.1
MT044257
MN996532
MG772934.1
MG772933.1
MT044258
MN996532
MG772934.1
MG772933.1
MT049951
MN996532
MG772934.1
MG772933.1
As suggested by the Z-curve visualization, three nearest SARS-like viruses are the same for all COVID-19 viruses.
Delete duplicates from the previous result and return their GenBank titles:
Our analysis shows that among 40 SARS-like viruses that we study here, Bat coronavirus isolate RaTG13, Bat SARS-like coronavirus isolate bat-SL-CoVZXC21, and Bat SARS-like coronavirus isolate bat-SL-CoVZC45 have closest genome similarity to COVID-19 viruses, providing evidence for the Bat origin of COVID-19.
Let us do some hierarchical clustering, too.
Concatenate Z-curves of COVID-19 and SARS-like viruses:
Create a numeric array of the type double-precision real our of each Z-curve:
Cluster the viruses based on their Z-curves using single-linkage clustering algorithm:
Assign different color to each cluster, and represent 3D images of corresponding Z-curves:
Project clusters into 3D using t-distributed stochastic neighbor embedding algorithm:
Constructs a dendrogram from the hierarchical clustering of the Z-curves (note the coloring is the same as the one for clusters and those that are bold and italic, are SARS-like viruses):

Concluding remarks

I generated Z-curves of 39 available genetic sequence of COVID-19 viruses together with 40 SARS-like viruses. An eyeballing of 3D visualization of Z-curves suggests that there are three SARS-like viruses very close to the virus family of COVID-19. By qualifying the similarity/dissimilarity between Z-curves, I identify those three viruses as Bat coronavirus isolate RaTG13 (MN996532), Bat SARS-like coronavirus isolate bat-SL-CoVZXC21 (MG772934 .1), and Bat SARS-like coronavirus isolate bat-SL-CoVZC45 (MG772933 .1). The hierarchical clustering of viruses based on the Z-curves also confirms the visualization, and identify Bat coronavirus isolate RaTG13 as the most-likely culprit of COVID-19. My results strongly support the hypothesis of a Bat origin of COVID-19.

References

R. Zhang and C-T Zhang, A Brief Review: The Z-curve Theory and its Application in Genome Analysis, Curr Genomics. 2014 Apr; 15(2): 78–94.
W-X Zheng et al, Coronavirus phylogeny based on a geometric approach, Molecular Phylogenetics and Evolution 36 (2005) 224–232