With more than 6.5 million infected people, SARS-CoV-2 is likely to have found a reservoir in the human population, increasing the possibility of an outbreak in the next autumn. To prevent the viral spread, an effective vaccine is needed.
More than 100 different vaccines are in development all over the word, including China, India, the USA and Europe. Most, if not all of them, target the spike protein using different protein formulations or platforms, such as a vectored vaccine or nucleic acids, RNA and DNA. The spike protein is the best candidate for historical reasons [1] , and for recent preclinical evidence in nonhuman primates [2] .
One of the main questions to ask is the following: from preclinical vaccine testing to the first-inman trial to confirmatory trials, regulatory approval and large-scale vaccine distribution, from 6 to 12 months will elapse. During this time frame, will the spike protein of the circulating virus mutate?
If so, are these genetic variants going to escape protective immune responses induced by the vaccine? In this paper, we try to provide a bioinformatic toolbox to address this potential issue.
A few preliminary reports already described the variant landscape of COVID-19 in different clades [3] [4] . Both pointed out the limited presence of functional mutations in key regions of the genome.
Nevertheless, some other preprints stressed a possible selection pressure on the D614G variant, the only frequent variation of the spike protein [5] . A recent bioinformatic analysis of SARS-CoV-2 epitopes showed a high homology with SARS-CoV [6] and, thus, characterized many possible Band T-cell epitopes, and the same group provided an in vivo characterization of the virus proteins that are most targeted by T cells, confirming the spike protein to be the first region of interest [7] .
To provide an evidence-based approach, we looked at the public available sequences and provide some comments with respect to vaccines and immune therapy. We also present a freely available bioinformatic protocol in order to enable all the bioinformatic community to dissect hCov-19 genomic variability and its functional impact on DNA vaccines. This work presents a protocol able to automatically mine a list of profiles useful for vaccine design from a pool of viral sequences and a target protein or domain of interest: a consensus sequence, in terms of average or minimal identity and a list of mutations with their absolute frequencies. Furthermore, the tool provides insightful graphical representation of the mutation frequency focusing on the genetic position of certain crucial domain, such as the receptor-binding domain (RBD) of the spike protein.
The overall workflow is depicted in Figure 1A . The only input of the analysis pipeline is a multisequence FASTA file representing the genome repertoire of many viral assemblies, which is available on the EpiCov™ section of the GISAID portal [8] . Filters used to fetch the 24,379 sequences analyzed in this manuscript were Host=Human, Virus=hCov-19, flagged "complete", "high coverage" and "low coverage exclusion".
As a first step, sequences are aligned against the Wuhan strain (NC_045512.2) with Bowtie2 [9] in order to produce a SAM/BAM alignment file. The said file is distributed across N files, containing one sequence each. This step is required for the bcftools [10] pileup command that computes all possible alleles found in every position, which is always with respect to the Wuhan strain reference. Variants are then called and indel-normalized via the bcftools call command. At this point, we have a list of genetic variants in a VCF file describing how every sample/viral individual differs from the original strain, and the next step would be describing and annotating the said variants in terms of functional effect and penetration in the viral population, while producing a consensus sequence that represents the most frequent allele for each position. The consensus command from bcftools builds the said consensus from the filtered VCF file. In addition to the consensus for the viral complete genome, we build the consensus sequence separately for the spike protein and RBD. Moreover, a consensus of minimal identity is constructed, applying the N character at variant sites instead of the most frequent base. Such consensus of minimal identity, by showing which residues are conserved and which residues are variable, would help in the design of a vaccine.
Variant annotation is employed via the snpEff package [11] that embeds the NC_045512.2 genome assembly annotation in its standard package. A human-readable (and easily parsable) table is formatted thanks to SnpSift jar package [12] of snpEff. 2D visualizations of the variants on the viral genome are generated via the LolliPlot function of the trackVignette's R package.
The whole workflow and test sequence files have been included in a Docker image, enabling users to run the whole pipeline with just one command (e.g. docker run covid-miner sequences.fasta).
The whole workflow, with sample test data and Dockerfile is available on https://gitlab.com/bioinfo-ire-release/covid-miner All software first picked for this workflow are standard tools employed in an NGS human sequencing facility lab. For instance, bowtie2 is a software design to handle short sequences (up to ~1000 bp according to the authors [13] ). Execution of the full pipeline on thousands of sequences resulted in 1-day computing time with 4-core and 16 GB RAM allocated. Testing and reimplementation of the workflow through optimized software, such as NUCMER [14] , which, from preliminary tests, cuts the computing time down to minutes, is ongoing at the time of the writing (Fig. 1A) .
The annotated variants table of all the GISAID sequences (available on June 3, 2020), is available in Table 1 . Looking at the total number of sequences the most frequent Spike protein variation is confirmed to be D614G ( Figure 1B and C) . After that, there is a set of variants with a frequency of 0.3-0.7% that occurred 100-200-times in the viral population. None of these variants lie in the Spike protein RBD locus. When focusing on the sole RBD, from amino acids 319-541, 13 variants arise, all with a relative frequency less than 0.1% and 10-20 absolute occurrences.
These RBD variants are represented in further detail in Figure 1B and C. We reasoned whether these variants belong to a certain geographical region and normalized the mutation count by "mutations per thousand isolates", that is
as the absolute mutation count and ‫ܫ‬ the total isolates for that country.
This normalization smooths out the inter-country variability, but the resolution bias remains due to the high variability in isolate sequencing among countries, changing in order of magnitude from thousands to dozens. This bias will cause rare mutations to be hidden in countries with a few associated isolates, while the associations with highly frequent outliers will remain more robust. By inspecting RBD variants in 12 countries, a few events, such as Furthermore, we explored whether these RBD mutations were influencing the class 1 presentation of SARS-CoV-2 epitopes. From the top five mutations at 0.1% frequency, we performed MHC class 1 binding prediction with NetMHC4 [15] , for both the wild-type (Wuhan strain) and the mutated strain, choosing the most frequent HLAs in public databases [16] . The software generated predictions for all possible 9-mers resulting from the full RBD sequence, and we computed how binders where calculated for all the considered HLAs. Out of 34 predictions, only three epitopes had binding affinity in mutated or wild-type, and none of them had an increased binding affinity caused by the mutation. Instead, one of them had a decreased binding affinity, namely GAEGFNCYF epitope caused by V483A mutation (Table 3) .
Finally, we asked whether these mutations had a significant impact on T-and B-cell recognition, and we ranked them via the Immune Epitope Database and analysis resource (IEDB) [17] . When
considering T-cells with the class 1 immunogenicity tool [18] , 117/250 (47%) epitopes showed a negative, non-immunogenic score in both wild-type and mutated forms. When focusing on the predicted class 1 binders for the same epitope/HLA combination in mutated peptides, both IYQAGSIPC and IYQAGNTPC epitopes do not show a strong induction in IEDB score, also starting with a low immunogenicity scores (Table 4 , Figure 2C ). These results point out to a limited variability in T-cell response mediated by these mutations.
Taken together, all these results confirm the both a limited variation pressure on critical SARS-CoV-2 loci, such as the RBD and on the whole Spike protein, and the scarce predicted immunological impact of the rare described variants.
This work presents a confirmatory study of SARS-CoV-2 genetic drift and a framework to dissect viral genomic variability by focusing on single mutations.
As a limitation, we acknowledge that the consensus sequence generated by our workflow do not represent any particular clade nor viral isolate and do not take into account linkage among variations. However, the focus in specific mutational events can enable an easier constant for a virus that is undergoing millions of replications, for clinical severity and vaccine efficacy monitoring.
Regarding the immunogenomic analysis, it would be of great interest coupling the class 1 prediction with B-rather than T-cell recognition. This task can be carried out via prediction tools, such as BepiPred [19] and is currently being investigated.
At the time of the writing, the clinical impact of the pandemic in Italy has been decreasing in the whole country in terms of contagion rates and mortality absolute counts. Several doctors working in COVID centers reported an anecdotical decreased severity of disease symptoms. Many of them explained these findings as a critical consequence of the change in clinical protocols, contact tracing, and higher temperatures that should decrease viral load and symptoms. In other cases, a molecular impact of mutations is being appointed as a warm gun. Our analysis, according to other reports, shows no substantial genetic pressure on the 83 available Italian genomes. 
Not applicable.
Not applicable.
All R and Bash/shell scripts developed are available at https://gitlab.com/bioinfo-ire-release/covidminer along with test files and a walkthrough.
The authors declare that they have no competing interests.
This work was supported by the Italian Ministry of Health (Ricerca Corrente 2019). of Country-RBD variant association C: Heatmap representing immunogenicity scores for wild-type and mutated RBD epitopes. 