CoVigator – a Knowledge Base for Navigating SARS-CoV-2 Genomic Variants

1. Introduction

The identification, characterization and monitoring of the pathogen responsible for a novel emerging disease is crucial for the development of a timely public health response. This includes rapid and open sharing of data [1] which has been adapted in past outbreaks to advance research and improve medical support [2,3]. The outbreak of the respiratory disease COVID-19 caused by SARS-CoV-2 demonstrated the increasing value of high-throughput sequencing by enabling the publication of the complete virus genome within one month of sampling [4,5]. The identification of the SARS-CoV-2 spike protein as a valuable target for vaccine design [6,7] led to the development of vaccines at unprecedented speed [8,9] and is still fostering further developments.

Nevertheless, the discovery of SARS-CoV-2 spike-glycoprotein mutants, associated with the potential to escape vaccine-induced protection, demonstrates the importance of monitoring SARS-CoV-2 genomic sequences to enable early detection of these genomic variants of concern. In a first study, we analyzed 1,036,030 genomic assemblies from the Global Initiative on Sharing Avian Influenza Data (GISAID) [10,11,12] and 30,806 Next Generation Sequencing (NGS) datasets from the European Nucleotide Archive (ENA). We reported non-synonymous spike protein mutations and their frequencies and analyzed the effect on known T-cell epitopes [13] . Although we confirmed low mutation rates of the spike protein, we experienced an increase in the number of genomic variants over time. Therefore, we further developed our pipeline to tackle potential escape mutations [14].

There are multiple initiatives to monitor SARS-CoV-2 mutations based on the GISAID dataset: NextStrain [15], CoV-GLUE [16] and Coronapp [17]. A further initiative based on the ENA dataset is the Galaxy project COVID-19 [18] ; other systems use regional data: CLIMB-COVID (COG-UK) [19] and CovRadar [20]. The COVID-19 Data Portal [21] is provided by EMBL-EBI and the European COVID-19 Data Platform to facilitate data sharing and accelerate research by making all the data available in the public domain and encouraging the research community to share SARS-CoV-2 data. Furthermore, there are some open source pipelines to identify mutations on SARS-CoV-2 data implemented in Nextflow: Cecret [22], nf-core viralrecon [23] and ncov2019-artic-nf [24].

Each dataset has its own advantages. While genomic assemblies are easier to share and interpret, raw reads provide granular information about the mutations through access to the pileup of reads supporting each mutation, also allowing the characterization of intrahost mutations. Analyzing both datasets together may be support the identification of potential false positives in the data and the confirmation of trends.

To enable monitoring of SARS-CoV-2 sequences from both sources, we have developed CoVigator, an NGS pipeline and dashboard that allows geographical and temporal navigation through SARS-CoV-2 genomic variants. We automatically download and analyze genomic assemblies from the COVID-19 Data Portal and raw reads from the European Nucleotide Archive (ENA). Furthermore, we screen the literature for studies on SARS-CoV-2 intrahost mutations [25,26] [27,28,29,30,31,32,33,34,35,36,37,38,39,40] and propose a filtering strategy to obtain a high quality set of intrahost mutations in the large and heterogeneous dataset obtained from ENA. Thus, the CoVigator platform supports early detection of variants, which potentially provides guidance for the further adjustment of vaccination strategies or therapeutics.

2. Materials and Methods

The CoVigator knowledge base is implemented in Python version 3.8 and the database for storing data is PostgreSQL (version 13.4).

The CoVigator pipeline (version 0.14.0) is implemented in the Nextflow framework version 19.10.0. All dependencies are managed within conda (version 4.9) environments [41] (see Table 1).

The CoVigator dashboard is also implemented in Python using the visualization framework Dash (version 2.1.0). The computation is distributed through a High Performance Computing cluster with a library that provides advanced parallelism, Dask (version 2022.9.2).

3. Results and discussion

System description

The CoVigator system (Figure 1) has three main components: 1) the knowledge base, 2) the analysis pipeline and 3) the dashboard. The knowledge base orchestrates for every sample the metadata retrieval, raw data download and finally its analysis through the pipeline for the detection of mutations. Furthermore it makes all necessary data available through a database (Postgre-SQL version 13). Finally, the dashboard presents the data to the end user through a set of interactive visualizations.

CoVigator operates via interaction with external systems: a high performance computing (HPC) cluster and the ENA and COVID-19 Data Portal Application Programming Interfaces (APIs).Samples between both original datasets (raw reads and genomic assemblies) may overlap. As recommended, some data providers might automatically upload both data formats. The results presented in the dashboard are stratified by dataset.

Knowledge base

The CoVigator knowledge base collects data from both genomic assemblies and raw reads, orchestrates its processing through the variant calling pipeline, and stores all the metadata, raw data and processed results in a relational database.

The data for both datatypes is fetched via the corresponding API hosted by the European Bioinformatics Institute [55], the metadata is normalized, the FASTQ (raw NGS reads) and FASTA (genomic assemblies) files downloaded and their MD5 checksums are confirmed to ensure data integrity.

Furthermore, the knowledge base iteratively builds a variant co-occurrence matrix (only for the raw reads dataset) and precomputes analyses on the data (binned abundance of mutations, dN/dS ratios per gene and domain, top occurring variants, pairwise co-ocurrence and counts of variants per lineage, country, sample, mutation type, length and nucleotide substitution) that ensure low latency responses.

Analysis pipeline

In general, the CoVigator pipeline processes FASTQ and FASTA files into annotated and normalized analysis-ready VCF files by two independent workflows (Figure 1 & Figure S1). We implemented the pipeline in the Nextflow framework [56] and managed all dependencies with Conda environments to enable seamless installation. We have embedded the SARS-CoV-2 reference genome ASM985889v3 [5]. Using a different reference, this pipeline could instantly analyze other virus sequences as well.

Dashboard

The dashboard is the user interface to CoVigator. There are two separate views for the raw reads and genomic assemblies datasets. Each view provides a set of tabs that allows the user to explore different aspects of the data held in the database. Each tab provides some interactive visualizations, described below. When applicable, the tabs provide a set of filters on the left side. These have been excluded from the screenshots for the purpose of clarity.

The most relevant tabs are described below and some notable findings are highlighted. The data shown here includes 137,025 samples downloaded from ENA on the 21^st October 2022; and 6,165,681 samples downloaded from the COVID-19 Data Portal on the 18^th November 2022.

Samples

The samples tab (Figure 2) enables the user to explore the accumulation of samples through time and the evolution of the dN/dS ratio on different genomic regions.

Figure 2A shows the accumulation of samples in each country. The dashboard allows the user to select specific countries and/or lineages.

Figure 2B,C shows the dynamics of the dN/dS ratio through time, over genes and protein domains. The dN/dS ratio aims to estimate the evolutionary pressure on the SARS-CoV-2 proteins and domains. This metric, although originally developed for assessing diverging species, is an imperfect but simple estimation of the evolutionary pressure within the same species [57,58], in this case SARS-CoV-2. There have been recent efforts to develop better alternatives for estimating the evolutionary pressure on SARS-CoV-2 [59]. The traditional interpretation of dN/dS is as follows: dN/dS < 0 indicates purifying selection, dN/dS = 1 indicates neutral evolution and dN/dS > 1 indicates positive selection.

Lineages

The lineages tab enables the user to explore the different lineages through time and geography (Figure 3). Both the accumulation of samples in every lineage worldwide (Figure 3A), and the dominant lineage through time (Figure 3B) can be viewed. In the screenshot, the displacement of B.1 by Alpha (B.1.1.7); then subsequently displaced by the multiple Delta lineages (AY.*); and finally displaced by the three Omicron lineages (BA.1, BA.2 and BA.3) can be seen.

Mutation statistics

The mutation statistics tab provides insights on the variant calling results on the different datasets and genomic regions (Figure 4). Expected trends in the data can be confirmed in these visualizations.

The median number of SNVs per sample in the raw reads dataset is 32, with an interquartile range (IQR) of 30 (Figure 4A). Additionally, the median number of MNVs is two with an IQR of one. The number of deletions is lower (median: 3, IQR: 2) than the number of SNVs and the number of insertions is even lower with few samples having just one insertion. For the genomic assemblies the numbers are slightly different with median SNVs 44 (IQR: 19), MNVs 2 (IQR: 0), deletions 4 (IQR: 3) and again just one insertion in few samples (Figure 4B).

We observe that the base substitution C>T is by large the most frequent, followed by G>T and A>G; the deletion TA>T and the MNV GG>AA is the most frequent in both datasets (Figure 4C,D).

In Figure 4E,F, we confirm that deletions are more frequent than insertions with an insertion-to-deletion ratio of 0.002 and 0.032 for raw reads and genomic assemblies, respectively. We also confirm two previous findings: 1) shorter deletions and insertions are more common than longer ones [60,61] and 2) the deletions and insertions not causing a frameshift are overrepresented as their impact in the resulting protein is more subtle [62]. In the genomic assemblies, we observe a long tail of deletions longer than 8 bp which is not observed in the raw reads results. We suspect this is a technical artifact introduced by our variant calling method. Finally, as shown in Figure 4G,H we observe that the most frequent mutation effect is a missense variant, followed by a synonymous variant. This is coherent between both datasets.

Recurrent mutations

The recurrent mutations tab allows the user to explore the most recurrent mutations by total count of observations through time within their genomic context (Figure 5).

In Figure 5A, the top recurrent mutations and their frequency and counts through time are shown. The size of the table can be parametrized for up to 100 mutations. For instance, the user can explore the most recurrent mutations in the whole genome, a given gene or a given protein domain. Furthermore, the period in which the monthly counts are shown can be parametrized. The gene viewer (Figure 5B) has multiple tracks: i) a scatter plot with the relevant mutations and their frequencies in the virus population, ii) ConsHMM conservation tracks and iii) gene and Pfam protein domains. The table in Figure S5 provides the decline and rise of the Alpha and Delta lineages, respectively, in the counts of mutations between April and July 2021.

Additionally, the mutation statistics tab provides a co-occurrence analysis that points to clusters of co-occurring mutations and their correspondence with virus lineages; or in the case of mutations shared between lineages, these clusters may contain a mixture of different but related lineages. Due to performance limitations, this analysis is only available in the raw reads dataset and at the gene level. In Figure S6, we show the Jaccard index co-occurrence matrix in the spike protein and its clustering results annotated with SARS-CoV-2 lineages in Table S1.

Clonal and intrahost mutations in the raw reads dataset

The FASTQ files provide the pile up of reads across the genome and this gives detailed information into the called variants. In particular, we can count the number of reads supporting each variant and this allows us to identify subclonal variants supported by only a fraction of the reads. These variants likely emerged within the host and are referred to as intrahost variants. The identification of intrahost variants is not possible on the genomic assemblies.

We consider high quality clonal mutations those with a VAF greater than or equal to 80 %, and those with a VAF greater than or equal to 50 % and lower than 80 % as low confidence clonal mutations. Only the high confidence clonal mutations are used to determine a consensus sequence and assign a SARS-CoV-2 lineage (Figure 6).

The remaining dataset of mutations poses a different technical challenge due to the difficulty to separate true low VAF mutations from noise. We first determine those mutations with a VAF below 50 % as raw candidate intrahost mutations.

We observed a large number of low frequency mutations among SARS-CoV-2 genomes. In order to establish a high quality set of intrahost mutations for studying viral evolution, we screened and compared the literature on SARS-CoV-2 intrahost mutations for different filtering approaches and implemented a conservative approach (Table 2).

4. Conclusions

The persistently increasing amount of publicly available SARS-CoV-2 sequencing data calls for robust platforms that allow constant monitoring of genomic SARS-CoV-2 variants in heterogeneous data sets. Our CoVigator pipeline covers the essential steps of preparing the data and calling variants from SARS-CoV-2 raw sequencing data from ENA and genome assemblies from the COVID-19 Data Portal. The pipeline is integrated within the CoVigator knowledge base that orchestrates download, processing and storage of the underlying samples and results. The CoVigator dashboard provides different visualizations and features for selecting clonal variants across all genes from the SARS-CoV-2 genome in a selected period. The dashboard also provides a comprehensive analysis of intrahost variants observed across detected mutations in the raw reads dataset. To this end, we propose a conservative filtering approach based on filtering samples and mutations. The dataset of intrahost mutations derived from public data that we make available through CoVigator is to the best of our knowledge the largest published dataset of SARS-CoV-2 intrahost mutations.

The identification of mutations over such heterogeneous datasets obtained with different sequencing protocols is challenging. With CoVigator, we observed VAF dilution on mutations identified by targeted amplicon sequencing with overlapping primers, genome edge effects and read edge effects. We aim to address these challenges in the future, e.g. by inferring the primers used in an arbitrary sample. Besides, we implemented a simplistic phasing of clonal mutations occurring in the same amino acid to ensure their correct annotation. However, we identified the need for a phasing method for low VAF mutations that existing germline phasing tools do not cover.

Future versions of CoVigator can be broadened to other use cases, such as other infectious organisms or co-existing infections during pandemic. Additionally, we envision the annotation of detected mutations to facilitate the selection of variants of concern for the development and evolution of preventive vaccines or therapies.