An efficient numerical representation of genome sequence: natural vector with covariance component

Natural vector

Natural vector is a 12-dimensional numerical representation of a nucleotide sequence (Deng et al., 2011), which is defined as follows. Let $S=s_1,s_2,s_3,\dots s_n$ be a genomic sequence of length n, and $L=\left\{A,C,G,\mathrm{T}\mathrm{o}\mathrm{r}\mathrm{U}\right\}$. For $k\in L$, the indicator functions $w_k(\cdot ):L\to \left\{0,1\right\}$ is defined as:

(1) $$\begin{array}{c}{w}_k\left(s_i\right)=\{\begin{array}{c}1,if{s}_i=k,\\ 0,otherwise.\end{array}\end{array}$$where $s_i\in L,i=1,2,3,\dots ,n.$ Let $n_k={\sum }_{i=1}^n{w}_k\left(s_i\right)$ denote the counts of nucleotide k in S, ${\mu }_k={\sum }_{i=1}^{n}i\frac{w_k\left(s_i\right)}{n_k}$ specify the average location of letter k, $D_j^k={\sum }_{i=1}^n\frac{{\left(i-{\mu }_k\right)}^j{w}_k\left(s_i\right)}{n_k^{j-1}{n}^{j-1}}$ be the j-th central moment of position of letter k. If $j=2$, we get the traditional 12-dimensional natural vector:

(2) $$\begin{array}{c}\left(n_A,n_C,n_G,n_T,{\mu }_A,{\mu }_C,{\mu }_G,{\mu }_T,D_2^A,D_2^C,D_2^G,D_2^T\right)\end{array}$$

Here we give an example to calculate the vector. If the genomic sequence is ACGGTAGTCC, the indicator functions are $w_A=1000010000,w_C=0100000011,w_G=0011001000,w_T=0000100100$. Each component of the vector is calculated as follow:

• • ${\mathrm{n}}_{\mathrm{A}}=2,{\mathrm{n}}_{\mathrm{C}}=3,{\mathrm{n}}_{\mathrm{G}}=3,{\mathrm{n}}_{\mathrm{T}}=2,$

• • ${\mu }_{\mathrm{A}}=1\cdot \frac{1}{2}+6\cdot \frac{1}{2}=3.5,$

• • ${\mu }_C=2\cdot \frac{1}{3}+9\cdot \frac{1}{3}+10\cdot \frac{1}{3}=7,$

• • ${\mu }_{\mathrm{G}}=3\cdot \frac{1}{3}+4\cdot \frac{1}{3}+7\cdot \frac{1}{3}=4.67,$

• • ${\mu }_T=5\cdot \frac{1}{2}+8\cdot \frac{1}{2}=6.5,$

• • $D_2^A=\frac{{\left(1-\frac{7}{2}\right)}^2}{2\cdot 10}+\frac{{\left(6-\frac{7}{2}\right)}^2}{2\cdot 10}=0.625,$

• • $D_2^C=\frac{{\left(2-7\right)}^2}{3\cdot 10}+\frac{{\left(9-7\right)}^2}{3\cdot 10}+\frac{{\left(10-7\right)}^2}{3\cdot 10}=1.27,$

• • $D_2^G=\frac{{\left(3-\frac{14}{3}\right)}^2}{3\cdot 10}+\frac{{\left(4-\frac{14}{3}\right)}^2}{3\cdot 10}+\frac{{\left(7-\frac{14}{3}\right)}^2}{3\cdot 10}=0.289,$

• • $D_2^T=\frac{{\left(5-\frac{13}{2}\right)}^2}{2\cdot 10}+\frac{{\left(8-\frac{13}{2}\right)}^2}{2\cdot 10}=\mathrm{0.225.}$

Then the 12-dimensional natural vector in formula (2) is:

$$\left(2,3,3,2,3.5,7,4.67,6.5,0.625,1.27,0.289,0.225\right).$$

A novel definition of covariance between nucleotides

The 12-dimensional natural vector definition reveals the respective distribution of four nucleotides, but ignores the relationship of each two nucleotides. We improve the 12-dimensional vector by defining the covariance between nucleotides. For the same genomic sequence $S=s_1{s}_2{s}_3\dots s_n$, and $L=\left\{A,C,G,\mathrm{T}\mathrm{o}\mathrm{r}\mathrm{U}\right\}$, we redefine the indicator function, $w_{kl}(\cdot )=w_{lk}(\cdot ):L\to \left\{0,1\right\},k,l\in L:$

(3) $$\begin{array}{c}{w}_{kl}\left(s_i\right)=w_{lk}\left(s_i\right)=\{\begin{array}{c}1,if{s}_i=korl,\\ 0,otherwise.\end{array}\end{array}$$

And the covariance between k and l is defined as:

(4) $$\begin{array}{c}Cov\left(k,l\right)=\sum _{i=1}^n\frac{\left[i-{\mu }_k\right]\left[i-{\mu }_l\right]w_{kl}\left(s_i\right)}{n\sqrt{n_k}\sqrt{n_l}}.\end{array}$$

The above formula reflects the correlation relationship of the position of each two nucleotides. The beauty of the definition is $Cov\left(k,k\right)=D_2^k$.

Thus, we get a novel 18-dimensional natural vector with covariance component:

(5) $$\begin{array}{c}\left(n_A,n_C,n_G,n_T,{\mu }_A,{\mu }_C,{\mu }_G,{\mu }_T,D_2^A,D_2^C,D_2^G,D_2^T,Cov\left(A,C\right),Cov\left(A,G\right),Cov\left(A,T\right),Cov\left(C,G\right),Cov\left(C,T\right),Cov\left(G,T\right)\right).\end{array}$$

For the same sequence ACGGTAGTCC, the indicator functions are $w_{AC}=1100010011,{\mathrm{w}}_{\mathrm{A}\mathrm{G}}=1011011000,w_{AT}=1000110100,w_{CG}=0111001011,w_{CT}=0100100111,w_{GT}=0011101100.$ The corresponding components are calculated as follows:

• • $Cov\left(A,C\right)={\sum }_{i\in \left\{1,2,6,9,10\right\}}^{}\frac{\left[i-3.5\right]\left[i-7\right]}{10\cdot \sqrt{2}\cdot \sqrt{3}}=2.06,$

• • $Cov\left(A,G\right)={\sum }_{i\in \left\{1,3,4,6,7\right\}}^{}\frac{\left[i-3.5\right]\left[i-4.67\right]}{10\cdot \sqrt{2}\cdot \sqrt{3}}=0.864,$

• • $Cov\left(A,T\right)={\sum }_{i\in \left\{1,5,6,8\right\}}^{}\frac{\left[i-3.5\right]\left[i-6.5\right]}{10\cdot \sqrt{2}\cdot \sqrt{2}}=0.85,$

• • $Cov\left(C,G\right)={\sum }_{i\in \left\{2,3,4,7,9,10\right\}}^{}\frac{\left[i-7\right]\left[i-4.67\right]}{10\cdot \sqrt{3}\cdot \sqrt{3}}=1.56,$

• • $Cov\left(C,T\right)={\sum }_{i\in \left\{2,5,8,9,10\right\}}^{}\frac{\left[i-7\right]\left[i-6.5\right]}{10\cdot \sqrt{3}\cdot \sqrt{2}}=1.735,$

• • $Cov\left(G,T\right)={\sum }_{i\in \left\{3,4,5,7,8\right\}}^{}\frac{\left[i-4.67\right]\left[i-6.5\right]}{10\cdot \sqrt{3}\cdot \sqrt{2}}=\mathrm{0.54.}$

Then the nucleotide distribution of sequence ACGGTAGTCC can be described by an 18-dimensional natural vector with covariance component:

$$\left(2,3,3,2,3.5,7,4.67,6.5,0.625,1.27,0.289,0.225,2.06,0.864,0.85,1.56,1.735,0.54\right).$$

The biological similarity between two sequences can be measured using the Euclidean distance of their corresponding 18-dimensional natural vectors, which is commonly used in our previous studies (Deng et al., 2011; Zhao, Tian & Yau, 2018; Zhao et al., 2019; Sun et al., 2021; He et al., 2020; Pei et al., 2020).

Datasets and tools

The newly proposed alignment-free method is tested on six microorganism datasets, including the genome and gene sequences of archaea, bacteria, virus and fungi, as shown in Table 1. The genomes of archaea, bacteria, virus and fungi are downloaded from National Center for Biotechnology Information (NCBI). Giant viruses in this study are regarded as the viruses with genome size greater than 170 kb, which are collected to study their relationship with archaea, bacteria, and other viruses. This giant virus dataset is collected from two public databases (Giant virus toplist: https://pitgroup.org/giant-virus-toplist/; The largest known viral genomes (completely sequenced, >170 kb): http://www.giantvirus.org/2014-04-14top.html) and virus dataset (Viruses with genome size greater than 170 kb). In order to compare with the previous nucleotide covariance study (Zhao, Tian & Yau, 2018), we re-download DNA barcode of fungi from Barcode of Life Data System (BOLD) according to the dataset record (Zhao, Tian & Yau, 2018).

To ensure the reliability of the data, we remove three types of sequences from the original datasets: (1) sequences without family taxonomic information; (2) families with less three sequences. (3) sequences of plasmid for bacteria dataset. DNA barcode uniquely identifies species using a short fragment of DNA sequence from specific genes (iBOL, 2022). For this DNA barcode of fungi dataset, we only remain the sequences pertaining to internal transcribed spacer (ITS) region of fungi (Conrad, Keith & Sabine, 2012; Naturvetenskapliga, 2010). At last, there are 298 genomic sequences belonging to 20 families for archaea dataset, 16,375 genomic sequences belonging to 178 families for bacteria dataset, 7,382 genomic sequences belonging to 83 families for virus dataset, 387 genomic sequences belonging to 22 families for fungi genome dataset, 677 sequences belonging to 16 families for giant virus dataset, and 95,542 sequences belonging to 467 families for fungi DNA barcode dataset. All accession numbers can be found in Datas S1–S6.

All the programs in this article are written in MATLAB R2020a and run on the same laptop (MacBook Air, 1.8 GHz Intel Core i5, 8 GB 1,600 MHz DDR3).

Convex hull principle for genomes

Convex hull principle for genomes has been demonstrated to be a wonderful classification tool (Zhao, Tian & Yau, 2018; Zhao et al., 2019; Sun et al., 2021). In Mathematics, the convex hull of a point set $A=\left\{a_1,a_2,\dots ,a_m\right\},a_i\in R^k$ is the minimal convex set that contains these m points, where $R^k$ is the k-dimensional Euclidean space, and $a_i$ is a k-dimensional point. The convex hull of finite set A is defined as the set of convex combinations of all points in A:

(6) $$\begin{array}{cc}CovA=& \{{\lambda }_1{a}_1+{\lambda }_2{a}_2+\dots +{\lambda }_m{a}_m:a_i\in A,{\lambda }_1+{\lambda }_2+\dots +\\ & {\lambda }_m=1,{\lambda }_i\ge 0,i=1,2,\dots ,m\}\end{array}$$

The boundary of the convex hull is spanned by some points of A (vertexes), and the rest points of A are lying inside the hull. The convex hull is a convex polygon if all $a_i\left(i=1,2,\dots ,m\right)$ are two-dimensional vectors, and the convex hull is a convex polytope if all $a_i\left(i=1,2,\dots ,m\right)$ are high dimensional vectors. Particularly, triangles and tetrahedrons are convex hulls.

Here $a_i\left(i=1,2,\dots ,m\right)$ are 18-dimensional natural vectors, and they are divided into several classes (families). Intuitively, the distribution of nucleotides from the same family is similar, and the 18-dimensional points from the same families lie closely. The sequence is transformed into an 18-dimensional numerical vector first, and the points from the same family can form a convex hull. Convex hull principle states that convex hulls corresponding to different families are mutually disjoint (Zhao et al., 2019).

Linear programming (LP) can be used to check whether two convex hulls are disjoint (Sun et al., 2021). If $A=Cov\left\{a_1,a_2,\dots ,a_m\right\}$ and $B=Cov\left\{b_1,b_2,\dots ,b_n\right\}$ intersect, then the convex combination of these points satisfy the formula: ${\sum }_{i=1}^m{\lambda }_i{a}_i={\sum }_{j=1}^n{\beta }_j{b}_j$, where ${\sum }_{i=1}^m{\lambda }_i=1,{\sum }_{j=1}^n{\beta }_j=1$. $a_i$ or $b_j$ is an 18-dimensional natural vector. It means that the following linear programming problem has a feasible solution (That is, there are non-zero coefficients $\left\{{\lambda }_1,{\lambda }_2,\dots ,{\lambda }_m;{\beta }_1,{\beta }_2,\dots ,{\beta }_n\right\}$ such that the minimum value of the optimization problem is 0):

(7) $$\begin{array}{c}min0\\ s.t.{\sum }_{i=1}^m{\lambda }_i{a}_i={\sum }_{j=1}^n{\beta }_j{b}_j\\ {\sum }_{i=1}^m{\lambda }_i=1,{\lambda }_i\ge 0,i=1,2,\dots ,m\\ {\sum }_{j=1}^n{\beta }_j=1,{\beta }_j\ge 0,j=1,2,\dots ,n\end{array}$$

The above problem can be implemented through linprog function built in MATLAB.

Convex hull classification of archaea, bacteria, virus and fungi

The five datasets, archaea, bacteria, virus, fungi (genome), fungi (DNA barcode), are used for convex hull classification, and the necessity of adding extra six-dimensional nucleotide covariance component to 12-dimensional natural vector is verified. For each dataset, we first calculate the 12-dimensional or 18-dimensional natural vector to describe the distribution of each sequence, and construct convex hull for each family. Then the LP method is utilized to check whether the convex hulls of different families are disjoint. The comparison of classification performance between 12-dimensional natural vector and 18-dimensional natural vector is displayed in Table 2. For archaea dataset, there are 190 convex hull pairs, and 184 convex hull pairs are mutually disjoint in 12-dimensional space, while all pairs are disjoint in 18-dimensional space. For bacteria dataset, there are $C_{178}^2=15,753$ convex hull pairs, and 15,160 pairs are not intersected using 18-dimensional vector method. The non-intersection ratio in 18-dimensional space is 96.24%, which is larger than that in 12-dimensional space (92.46%). For virus dataset, there are 3,403 convex hull pairs. A total of 3,322 pairs are disjoint in 18-dimensional space, and 3,321 pairs are disjoint in 12-dimensional space. For both fungi (genome) and fungi (DNA barcode) datasets, the non-intersection ratios in 18-dimensional space are higher than those in 12-dimensional space. The results demonstrate that our new 18-dimensional natural vector performs better classification results than 12-dimensional natural vector. The 18-dimensional vector contains more sequence information because it considers the correlation relationships of the two nucleotides. The results of fungi (DNA barcode) dataset further show that covariance can be used not only for sequence classification at genome level, but also for sequence classification at gene sequence level.

To visualize the convex hull classification results, we use the traditional support vector machine method (SVM) to reduce the dimension (Sun et al., 2021). For two disjoint convex hulls in 18-dimensional space, there is a hyperplane to separate them:

(8) $$\begin{array}{c}{w}^{T}x+b=0,\end{array}$$ $w={\left(w_1,w_2,\dots ,w_k\right)}^T$ is the normal vector, b is the offset item. There is at least a vector $v$ perpendicular to $w$. For an 18-dimensional point, named $NV$ in the convex hull, the new point projected into two-dimensional space is $\left(w^T\cdot NV,v^T\cdot NV\right)$. Then the convex hull can be visualized in two-dimensional space. As shown in Fig. 1, we randomly select four disjoint convex hull pairs. Points from the same family gather together.

Preliminary statistical analysis of the genomes of giant viruses

The research on phylogenetic relationship between giant virus and bacteria, archaea and other viruses is still active (Claverie & Abergel, 2013). The discovery of Mimivirus comes as a shock in biological field for a long time (Birtles, Rowbotham & Storey, 1997), and its genome is about 118 kbp (Didier, Stéphane & Catherine, 2004). Later viral genomes larger than 118 kbp started to accumulate quickly. This suggests that viruses with genomes larger than 118 kbp are not as exceptional as previously thought, especially Pandoravirus, which is discovered in recent years, with a length of more than 2,500 kbp (Philippe et al., 2013). Recent discoveries are not sufficient to erase the tremendous gap in genome size between giant viruses (Brandes & Linial, 2019), which is also reflected in Fig. 2. The x-axis represents the sequence number, the y-axis represents the genome size, and each dot represents a sequence. There are many gaps between genomes sizes of “Giant Virus” (blue color). Here “Giant Virus” only indicates the sequences of giant virus collected from two public databases (genome size > 170 kb). “Giant & RefSeq” indicates the sequences of giant virus collected from virus dataset (genome size > 170 kb). “Other Virus RefSeq” indicates the rest sequences in virus dataset except for giant virus. In Fig. 2A, there is a significant slope change at about 300 kbp, and the genome size of giant viruses is larger than that of other viruses (green color). In Fig. 2B, it is found that the genome sizes of several giant viruses and bacteria overlap at about 2,000 kbp, and the genome sizes of some viruses are larger than those of some bacteria. In Fig. 2C, the genome sizes of several giant viruses also overlap at about 2,000 kbp. Fig. 2 gives an intuitive understanding of the genome size relationship between giant virus and bacteria, archaea, other viruses. There is not yet any reason to believe that the genome size of giant viruses has reached the upper limit, and viruses with different genome lengths may emerge continuously (Claverie & Abergel, 2013). The origin of the giant virus still remains mysterious, so we can study it through the relationship between genomes of similar length.

We next make preliminary statistics on the gene distribution to understand the giant virus. All known giant viruses belong to the phylum Nucleocytoviricota (nucleocytoplasmic large DNA viruses, NCLDVs), which contains many unique genes that cannot be found in other life forms (Ogata, Toyoda & Tomaru, 2009). Here we only consider the coding sequence. We plot the number of NCLDV coding sequences on the y-axis, and the genome size on the x-axis. The results are shown in Fig. 3. The representative families of giant viruses, Mimiviridae and Pandoraviridae have larger genome size and more coding sequences (Mimiviridae: purple star, Pandoraviridae: black triangle). In addition, the G+C content of Pandoraviridae is the highest among NCLDVs (average 61.3%), and the earliest discovered virus group, Mimiviridae has the lowest G+C content (average 26.3%, Table S1).

Nearest neighbor classification of the genomes of giant viruses

The similarity of genomic distribution between giant viruses and bacteria, archaea, other viruses is of great importance to understand their phylogenetic relationship. We use a new 18-dimensional natural vector with covariance to measure it. There are 24,250 different biological sequences in these four datasets (Genomes of giant virus, bacteria, archaea, and virus), of which 677 sequences belong to 31 groups for giant virus dataset. Each sequence is converted into an 18-dimensional natural vector first. For each sequence in giant virus dataset, the closest sequence is then calculated and the sequence type is checked. The Euclidean distance between two 18-dimensional natural vectors is used to measure the biological similarity of the corresponding sequences. The nearest neighbor classification results are shown in Table 3.

Except for Polydnaviridae and Myoviridae, the nucleotide distribution of most sequences is similar to that of bacteria $\frac{285}{309}\times 100\mathrm{\%}=92.23\mathrm{\%}$. About half of the 138 Myoviridae sequences are close to the bacterial genomes (67 sequences) and half to the viral genomes (71 sequences). One of the reasons for this phenomenon is that the sequence length of Myoviridae is similar to that of some bacteria. Of the 230 Polydnaviridae sequences, 229 are closest to the virus sequence. These 230 sequences belong to two genera, Ichnovirus (IV) and Bracovirus (BV), and five species (Table S2). Glypta fumiferanae ichnovirus (GfIV) has 105 segments, which is the virus with the most segments among all segmented viruses. The segment length is about 2,777 bp, which is shorter than bacterial sequence. In addition, all sequences of three representative giant virus families, Mimiviridae, Marseilleviridae, and Pandoraviridae are closer to bacterial sequences.

Phylogenetic analysis of the genomes of giant viruses

Phylogenetic relationships between giant virus and bacteria, archaea, other viruses are studied using our new alignment-free method. We select three representative families from giant virus (Mimiviridae, Pandoraviridae, and Marseilleviridae) and randomly select two or three families from other three datasets to construct the phylogenetic tree. Figure 4 shows the phylogenetic tree of three other virus families (Adenoviridae, Anelloviridae, and Closteroviridae) and the three families of giant virus, different groups are marked in different colors. All groups are separate except Mimiviridae. To clearly observe the relationships between the three giant virus families, two bacterial families, Pseudomonadaceae and Alteromonadaceae are selected for further analysis. The phylogenetic result in Fig. 5 shows that all sequences from the same family are clustered except Mimiviridae. Beyond that, two archaea families, Methanosarcinaceae and Natrialbaceae are selected to construct the phylogenetic tree, the same family groups together except Mimiviridae is shown in Fig. 6. The algorithm of constructing the phylogenetic tree is unweighted pair-group method with arithmetic mean (UPGMA), which is an approach of constructing rooted tree based on distance matrix. The distance matrix is calculated using the commonly used Euclidean distance. In Figs. 4–6, the same 10 Mimiviridae species and Marseilleviridae species are clustered together. The sequence information is displayed in Table S3. The length of the 10 sequences is shorter than that of other viruses of Mimiviridae, but similar to some viruses of Marseilleviridae. The average sequence length of Marseilleviridae is about 370 kbp, that of Mimiviridae is about 1,000 kbp, and the length of 10 sequences is hundreds of thousands of bp.

To demonstrate the importance of natural vector with six-dimensional covariance component in characterizing biological sequence, we compare it with the traditional 12-dimensional natural vector and another alignment-free method, position natural vector (He et al., 2020) on the same three small datasets. The phylogenetic trees based on 12-dimensional natural vector are shown in Figs. S1 to S3. In Fig. S1, the sequences of Adenoviridae colored navy blue are not clustered together. In Fig. S2, there are two extra sequences from Mimiviridae clustered with Marseilleviridae. In Fig. S3, the sequences of Methanosarcinaceae are not grouped together. Additionally, the results of position natural vector method are shown in Figs. S4 to S6. In Fig. S4, a sequence from Marseilleviridae is incorrectly clustered with Mimiviridae. Figure S5 shows an unstable result compared with Fig. 5. In Fig. S6, the sequences of three families Methanosarcinaceae, Pandoraviridae, and Natrialbaceae are mixed. In conclusion, the tree based on our new 18-dimensional natural vector gives more reasonable and interpretable results.

Haursdorff distance can reflect the genetic distance between groups, and we want to analyze the interspecific differences of the above 10 families. Mathematically, Haursdorff distance is used to calculate the distance between point sets. If X and Y are two non-empty point sets, then the Haursdorff distance is defined as:

(9) $$\begin{array}{c}d\left(X,Y\right)=max\left\{\underset{x\in X}{sup}\underset{y\in Y}{inf}d\left(x,y\right),\underset{y\in Y}{sup}\underset{x\in X}{inf}d\left(x,y\right)\right\}.\end{array}$$

The distance calculation can be implemented through the code of the MathWorks community: https://ww2.mathworks.cn/matlabcentral/fileexchange/26738-hausdorff-distance. Each sequence is transformed into an 18-dimensional natural vector, then Haursdorff distance between groups is calculated, and UPGMA algorithm based on distance matrix is used to construct the tree, as shown in Fig. S7. Mimiviridae is marked with circles, and Marsellieviridae with triangles. Mimiviridae is closer to the root of the tree than Marsellieviridae.

Time and phylogeny comparison with previous covariance method

The previous study has shown that nucleotide covariance is related to the phylogenetics of fungi (Zhao, Tian & Yau, 2018), we redefine the nucleotide covariance to make it simpler and more general.

The previous covariance definition is:

(10) $$\begin{array}{c}Cov\left(k,l\right)=\frac{Cov\left(A,B\right)}{N},\end{array}$$

$k,l\in \left\{A,C,G,\mathrm{T}\mathrm{o}\mathrm{r}\mathrm{U}\right\}.A=\left\{a_1,a_2,\dots ,a_n\right\},B=\left\{b_1,b_2,\dots ,b_m\right\}$ are the position of k and l, respectively.

• • If $m=n$, $Cov\left(A,B\right)={\sum }_{i=1}^n\frac{\left(a_i-{\sum }_{i=1}^n\frac{a_i}{n}\right)\left(b_i-{\sum }_{i=1}^m\frac{b_i}{m}\right)}{n}$.

• • If $m\ne n$ (assume $m>n$), the covariance between A and any n values in B is computed and the average of these $C_m^n$ results is $Cov\left(A,B\right)$.

For a bacterial sequence (Accession number in Genbank is AP018515, sequence length is 3,041,114 bp), the number of nucleotide T is m = 723,764 and the number of nucleotide C is n = 799,865, the covariance of T and C using the previous definition is:

(11) $$\begin{array}{c}Cov\left(T,C\right)=\frac{Cov\left(A,B\right)}{3041114},\end{array}$$

Suppose the positions of T and C are $A=\left\{a_1,a_2,\dots ,a_{723764}\right\},B=\left\{b_1,b_2,\dots ,b_{799865}\right\}$, respectively, then the covariance between A and any 723,764 values in B is computed and the average of these $C_{799865}^{723764}$ results is:

(12) $$\begin{array}{cc}Cov\left(A,B\right)=& \frac{1}{C_{799865}^{723764}}{\sum }_{\{j_1,j_2,\dots ,j_{723764}\}\in B=\{b_1,b_2,\dots ,b_{799865}\}}\\ & {\sum }_{i=1}^{723764}\frac{\left(a_i-{\sum }_{i=1}^{723764}\frac{a_i}{n}\right)\left(b_i-{\sum }_{i=1}^{723764}\frac{b_{j_i}}{m}\right)}{723764},\end{array}$$

If the number of two nucleotides is quite different, the above formula (12) is too complicated and difficult to compute, which will take a lot of time to calculate. Our new nucleotide covariance definition is $Cov\left(k,l\right)={\sum }_{i=1}^N\frac{\left[i-{\mu }_k\right]\left[i-{\mu }_l\right]w_{kl}\left(s_i\right)}{N\sqrt{n_k}\sqrt{n_l}}$, there is no need to compute $C_{799865}^{723764}$ values, so the calculation of covariance takes less time.

We compare the running time of the two definitions using a sequence segment with a length of 90 bp (gaggaagtaaaagtcgtaacaaggtttccttccgggtgtagcacctgccgaagcctcccgcagcgactctaaagaaactgcgcagtctgc, a segment of a sequence whose accession number is DQ525472). The results are displayed in Fig. 7. The calculation time increases greatly with the increase of the sequence length for the previous covariance method, but is robust for our new method. It costs about 0.5 s to compute the previous covariance vector of a sequence with length 75 bp, which is much longer than that of our method (Fig. 7). While it only costs 0.017 s to compute our covariance vector of a sequence with length 813 bp (GenBank accession number is DQ525472).

The previous covariance method is tested on the fungi DNA barcode dataset, so we also apply our method on the same dataset. To show that our method can also be used for phylogenetic analysis, we choose 311 nucleotide sequences of DNA barcode from nine species of fungi (Monosporascus cannonballus, Peniophorella praetermissa, Biatora helvola, Ceratobasidium cereale, Rhodotorula glutinis, Geosmithia landonii, Fusarium cf. solani, Acarospora smaragdula, Acaulospors kentinensis) to construct phylogenetic tree, as shown in Fig. S8. The nine groups can be separated based on our method, while a sequence from Peniophorella praetermissa doesn’t cluster together with Peniophorella praetermissa based on the previous covariance method.