Unsupervised learning analysis on the proteomes of Zika virus

Data preprocessing

The whole flow chart is shown in the Fig. 1. It started by retrieving a total of 391 nonredundant ZIKV proteomes along with their metadata (available as Supplementary Materials) for hosts (human, mosquito, and primate) and geographic regions from the virus variation database (https://www.ncbi.nlm.nih.gov/genomes/VirusVariation/Database/nph-select.cgi?cmd=database, accessed on 08-09-2023) (Brister et al., 2013). After eliminating proteomes with unambiguous amino acids, the final data set was composed of 293 proteomes. This sequence data was employed to construct multi alignments with the MAFFT (Multiple Alignment using Fast Fourier Transform) software (Katoh, Rozewicki & Yamada, 2019) using default options (available as Supplementary Materials). The proteome phylogenetic tree was constructed from the multi-alignment with IQ-TREE v. 2.3.6 software, which select the most appropriate amino acid evolution mathematical model (Minh et al., 2020) (http://www.iqtree.org/). The resulting newick format tree file was plotted and colored based on the geographic metadata with the iTOL web platform (https://itol.embl.de/) (Letunic & Bork, 2021) (available Supplementary Materials).

Description of used models

URF: The ML algorithms were implemented in house-scripts running on R 4.2.2 version (R Development Core Team, 2016) as follows: First the amino acid polymorphic sites were extracted from the multi-alignments with the function “alignment2genind” from the adegenet R package (Jombart, Balloux & Dray, 2010), from which 422 polymorphic sites were determined. Based on these data, a synthetic data set of the same size was produced by permutating the rows of each column and was then combined with the original data. This information was used to perform a Random Forest (RF) analysis with the R package randomForest using the function “randomForest” (Breiman, 2001; Liaw & Wiener, 2002). The hyperparameters were optimized using a five-fold cross-validation with three-repeat procedure with the caret R package (Kuhn, 2008). The “mtry” and “ntree” parameters were settled using grid search, that include the square root (sqrt) of the predictors and 500, 1,000, 1,500, and 2,000 trees. The RF out of bag (OOB) error was used as partial validation of goodness of fit for the URF. The final best hyperparameters were then used to perform the URF on the original data by including the parameter proximity = TRUE, which allowed to obtain the similarity (or proximity) matrix. The produced matrix was converted to a dissimilarity matrix by subtracting 1 for each data to compress all the information as previously described (Liaw & Wiener, 2002; Afanador et al., 2016). Then the matrix was converted to a Euclidean distance matrix to confirm the structure on the data via heatmap plotting and the Hopkins statistic. Finally, the dissimilarity matrix based on real data was further used as input to the PCA, t-SNE, UMAP, and autoencoders dimensionality reduction algorithms. The pseudocode that includes the URF and the DR procedures is shown below.

1. Combined_Data=real_data + synthetic_data,

2. Unsupervised RandomForest(data=combined_data, tree=num_trees, mtry=min_samples_split, proximity= proximity_matrix):

3. forest = list()

4. proximity = matrix(len(data), len(data))

5. For i in tree do:

6. bootstrap_sample = bootstrap(data)

7. tree = randomForest(bootstrap_sample, mtry)

8. forest.aggregate (tree)

9. For each pair of samples (x_i, x_j) in data do:

10. If x_i and x_j end up in the same leaf node in the tree then:

11. proximity_matrix[i][j] += 1

12. End If

13. End For

14. return (forest, proximity_matrix)

15. return (proximity_matrix/ntree)

16. evaluate error rate on OOB

17. UnsupervisedRandomForest(data=real_data, tree=num_trees, mtry=min_samples_split, proximity=TRUE):

18. distance_matrix = 1 - proximity_matrix

19. cluster with PCA(distance_matrix)

20. cluster with t-SNE(distance_matrix)

21. cluster with UMAP(distance_matrix)

22. cluster with Autoencoders(distance_matrix)

PCA: The PCA was applied using the function “prcomp” and the option “scale=FALSE” with the R base package. PCA uses scaled data to calculate the covariance matrix to identify correlated variables. Then this information is used to derive new uncorrelated variables called the principal components that retain the variance of the original data into few axes’ coordinates.

t-SNE: The t-SNE algorithm was applied with “Rtsne” function of the Rtsne package by settling the perplexity parameter at 50 with 1000 step iterations to identify the similar and dissimilar data points. t-SNE reduces the high dimensional data using a symmetrized cost function and Student-t distribution to calculate the similarity in the low-dimensional space of two data points randomizing the initialization (Van der Maaten, 2014; Krijthe, 2015). To balance local and global features in the data, t-SNE uses the perplexity parameter that varies in the range of 5 to 50 to guess the closest neighbors for each data point.

UMAP: UMAP was implemented with the “umap” function from the umap R package that considers 200 epochs and 15 neighbors as cut-off to compute the similarity. The UMAP algorithm is similar to the t-SNE and uses the t-distribution to calculate the similarity in the data except that UMAP initializes in the same data points through spectral embeddings that allows to move a single or group of data points at once to calculate its similarity, making it faster than t-SNE (McInnes et al., 2018).

AutoEncoders: Autoencoders neural network was constructed using the Keras package for R (Chollet, 2015). The neural architecture was designed containing three layers with 12 neurons in the first encoder layer, six in the second bottleneck layer and 12 in the third decoder layer using the activation hyperbolic tangent function settled with the parameter “tanh”. Using these parameters the AE shrinks the data into the bottleneck with 2,000 training iteration epochs.

Validation

External validation was carried visually by labeling the recovered clusters with the metadata information through two-dimensional plots resolved by each DR; thus, if the clusters formed according to its geography and host isolation, they were considered true clusters. In addition, the internal validation was performed through silhouette plots using k-means and the Euclidean distance in the elbow analysis to identify the number of suitable groups produced by each DR. The silhouette coefficient was considered the mean distance between samples in specific groups and other clusters. Thus, the larger the coefficient, the more consistency is in the cluster distribution. Finally, the whole unsupervised analysis strategy was applied to a larger dataset retrieved from virus variation (https://www.ncbi.nlm.nih.gov/genomes/VirusVariation/Database/nph-select.cgi?cmd=database, accessed on 26-07-2024) comprising 1,343 whole dengue type 2 (DENV2) proteomes isolated from human hosts (available as supplementary materials). As these DR methods could produce slightly different results in each training test, six test runs were performed and analyzed to confirm the reproducibility of the cluster patterns.

Phylogenetic correlation and dependence test analysis

The R packages ape and adephylo (Jombart, Balloux & Dray, 2010) were used to read and merge the axes cluster coordinates with the phylogeny using the function “phylo4d”. Then the phylogenetic correlations and statistical test were calculated with the package phylosignal (Keck et al., 2016). This package contains several functions to calculate correlations and statistical test of phylogenetic dependence of continuous trait values with the phylogenetic branches. The phylogenetic correlograms were produced with the function “phyloCorrelogram”. The Abouheif’s Cmean and Pagel’s Lambda tests were calculated with the functions “phyloSignal” specifying the methods “Cmean” and “Lambda” for each test, respectively. A p value <0.01 was deemed as significant in the phylogenetic statistical test under the null hypothesis of lack of phylogenetic dependence.

Unsupervised random forest retrieves cluster tendency

Using the best hyperparameters (mtry=10, ntree=1,500, Fig. S1) the RF indicates an OOB of 1.19%, confirming that it can distinguish the real data from its synthetic version. Moreover, the Hopkins statistic value was of 0.91, indicating the existence of clustering tendency in the real data. The visual exploration confirmed the cluster tendency which was according to the geographic origin (Fig. 2). In the heatmap six well-defined blocks are formed, using the Euclidean distance.

UMAP showed the best performance for cluster separation

To prove the clustering performance for PCA, t-SNE, UMAP and Autoencoders, the k-means clustering using the elbow curve analysis based on silhouette coefficient was applied. The analysis of the elbow curve showed similar separation distances (average silhouette coefficient >0.6, Figs. 3A, 3D, 3E and 3G) among the four tested DRs, but different number of clusters of k (range 3 to 8). Among them, PCA and UMAP were the best showing the ability to identify more clusters than the other algorithms (Figs. 3B, 3D, 3F and 3H), but UMAP showed better average silhouette coefficient and cluster separation than PCA (Figs. 3E and 3F).

Geography is the main factor for cluster separation

After labeling the clusters, the four algorithms showed that ZIKV isolates were grouped according to their host (Fig. 4). For example, those from arboreal mosquitoes clustered with NHP ZIKV isolates and the Aedes aegypti isolates clustered mainly with human isolates. In the cluster analysis, the African ZIKV group was clearly separated from the other geographical regions, as it consisted of NHP and arboreal mosquitoes and the human isolate (virus variation ID AMR68906), supporting the view of limited adaptation of sylvatic strains compared to endemic urban strains (Beaver et al., 2018). Interestingly, a clear geographical dispersion pattern was observed for some Asian ZIKV that tend to migrate into the European, Oceanic, North American, and South American regions. Furthermore, these Asian isolates tend to mix tightly with some American isolates. The t-SNE and UMAP and showed the clearest separation within regional clusters, in which the North American cluster harboring few Asian genomes was of note (Figs. 4D and 5F). In this cluster, the DR analyzes identified three Asian strains (sequence IDs ASU55505, ATI21641, and ARI68105), which were previously identified as imported American strains in Hainan China and Singapore (Tan et al., 2018).

The analytical strategy was generalizable to an external proteome dataset

To confirm whether the analytical strategy for unsupervised analysis takes a broad view to other real datasets, we applied the whole analytical strategy to a data matrix of polymorphic sites extracted from a multialignment of 1,343 DENV2 proteomes. We found that the dissimilarity matrix produced by the URF recovers the cluster tendency (Fig. S2) with a Hopkins statistic of 0.93. The UMAP was also the best DR algorithm compared to the others (Fig. S3, Fig. 5C). External validation by labeling with metadata also showed that geography is the main factor shaping the DENV2 dispersion (Figs. 5A, 5B, 5C and 5D), but was particularly evident in the UMAP (Fig. 5C). This result confirmed that this analytical strategy could be generalized to analyze the dispersal of viral pathogens.

Phylogeny correlates with cluster coordinates

As expected, the phylogeny showed similar clustering tendencies; for instance, the African ZIKV were clearly separated from the other groups with a larger branching length. The same Asian ZIKV were mixed with European, Oceanic, North American, and South American regions with smaller branching length (Fig. 6).

The phylogenetic clustering was related to the first two-dimension coordinates obtained by the four DRs. A phylogenetic correlogram analysis showed to some extent significant positive correlations of the coordinates with smaller branch lengths and significant negative correlations with larger branch length (Figs. 7A, 7B, 7C and 7D). Additionally, the Abouheif’s Cmean and Pagel’s Lambda tests (Table 1) showed also significant phylogenetic dependency with the dimension coordinates. These results confirm that UL retains the evolutionary findings as the phylogenetic gold standard methods without the need of an underlying evolutionary method.