SARS-CoV-2 CoCoPUTs: analyzing GISAID and NCBI data to obtain codon statistics, mutations, and free energy over a multiyear period

Abstract

A consistent area of interest since the beginning of the severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) pandemic has been the sequence composition of the virus and how it has changed over time. Many resources have been developed for the storage and analysis of SARS-CoV-2 data, such as GISAID (Global Initiative on Sharing All Influenza Data), NCBI, Nextstrain, and outbreak.info. However, relatively little has been done to compile codon usage data, codon-level mutation data, and secondary structure data into a single database. Here, we assemble the aforementioned data and many additional virus attributes in a new database entitled SARS-CoV-2 CoCoPUTs. We begin with an overview of the composition and overlap between two of the largest sources of SARS-CoV-2 sequence data: GISAID and NCBI Virus (GenBank). We then evaluate different types of sequence curation strategies to reduce the dataset of millions of sequences to only one sequence per Pango lineage variant. We then performed specific analyses on the coding sequences (CDSs), including calculating codon usage, codon pair usage, dinucleotides, junction dinucleotides, mutations, GC content, effective number of codons (ENCs), and effective number of codon pairs (ENCPs). We have also performed whole-genome secondary RNA structure prediction calculations for each variant, using the LinearPartition software and modified selective 2'-hydroxyl acylation analyzed by primer extension (SHAPE) data that are available online. Finally, we compiled all the data into our resource, SARS-CoV-2 CoCoPUTs, and paired many of the resulting statistics with variant proportion data over time in order to derive trends in viral evolution. Although the overall codon usage of SARS-CoV-2 did not change drastically, in line with the previous literature on this subject, we did observe that while overall GC% content decreased, GC% of the third position in the codon was more positive relative to overall GC% content between February 2021 and July 2023. Over the same interval, we noted that both synonymous and nonsynonymous mutations increased in number, with nonsynonymous mutations outpacing synonymous mutations at a rate of 3:1. We noted that the predicted whole-genome secondary structures nearly all contained the previously described virus-activated inhibitor of translation (VAIT) stem loops, validating for the first time their existence in a whole-genome secondary structure prediction for many SARS-CoV-2 variants (as opposed to previous local secondary structure predictions). We also separately produced a synonymous mutation-deprived set of SARS-CoV-2 variant sequences and repeated the secondary structure calculations on this set. This revealed an interesting trend of reduced ensemble free energy compared to the unaltered variant structures, indicating that synonymous mutations play a role in increasing the free energy of viral RNA molecules. These data both validate previous studies describing increases in viral free energy in human viruses over time and indicate a possible role for synonymous mutations in viral biology.

Keywords: SARS-CoV-2; VAIT; bioinformatics; codon usage; secondary structure.

Figure 1.

Beginning with raw sequence data and metadata from GISAID and NCBI, we QC sequences using Nextclade software, split sequences into individual FASTA files for each Pango variant, identify the overlap (matched metadata and sequence identity) between GISAID and NCBI to create a deduplicated combined sequence resource, use a pseudo-ancestral sequence curation strategy to reduce the large cohort of sequences to one sequence per variant, and finally perform our sequence composition and secondary structure/free energy calculations for presentation on our website and within this manuscript.

Figure 2.

(a) Five different categories of “matching” between NCBI and GISAID sequences prior to QC, ranging from a complete sequence and metadata match to a complete sequence and metadata mismatch; (b), the regional distribution (by continent) of the entire combined dataset; (c–d) the regional distribution of sequences from GISAID, NCBI, and the overlap (only) between the two.

Figure 3.

(a) The submissions per month from GISAID, NCBI, and the overlap between them; (b) The monthly submissions by country, where only the top 10 countries by total volume of submissions are included, while all other countries are grouped into the “All Others” category; (c) The total submissions by country, again showing only the top 10; (d) The same countries as a percentage of the total human population (based on publicly available UN data).

Figure 4.

(a) The initial landing page, with the option to search multiple SARS-CoV-2 strains in a single query or split them across multiple queries; the former will result in their results being aggregated together, whereas the latter will result in them being compared against each other; (b) One of the “Results” panels, which in this case shows a bar graph of codon usage between two different queries (e.g. lineage “A.2 and AL.1” and lineage “AY.4.”

Figure 5.

(a–c) Respectively, the percent abundance of each individual codon, dinucleotide, or junction dinucleotide was first calculated for each variant, after which the difference in abundance was calculated compared to the wild-type reference, and these difference values were computed for each time point based on the variant composition in that month; (d) The difference in overall GC content (GC%) as well as the difference in GC content for the first, second, and third positions within each codon (GC1–3%) compared to the wild-type for each time point, based on the variant proportions at that time point.

Figure 6.

(a-b) Changes relative to wild-type in CAI and CPAI (respectively) over a 46-month period for both the whole-genome CDS (first row) and all individual genes (subsequent rows); (c) Changes relative to wild-type in ENC, ENC with GC correction, ENCP, and ENCP with GC correction over the same time period; for all panels, red color is used to depict points where the CAI or CPAI exceeded the maximum range of the color bar (only applies to ORF8).

Figure 7.

(a) Changes in free energy over a 46-month period for both unmodified SARS-CoV-2 variants and a special “synonymous removed” set of variant sequences; (b) changes in both synonymous and nonsynonymous mutations over time; (c) autocorrelation statistic for the unmodified free energy time series.

Figure 8.

(a—f) The number of nucleotides attributable to the 5ʹ unpaired end, stems, interior loops, multiloops, hairpin loops, and the 3ʹ unpaired end (respectively) over time.

Figure 9.

(a) The secondary structure conservation in our cohort of thousands of SARS-CoV-2 variants in our dataset; (b) the number of stem-loops identified over time over a 46-month period, with each time point representing the number of stem-loops for the unique distribution of SARS-CoV-2 variants at that time.