VOE: automated analysis of variant epitopes of SARS-CoV-2 for the development of diagnostic tests or vaccines for COVID-19

Experimental setup

All experiments in this study were performed in a 64-bit architecture server with a 1 TB SSD, 64 GB of RAM, two Intel(R) Xeon(R) Silver 4110 CPUs @ 2.10 GHz running Ubuntu 22.04.2 LTS. All analysis software was installed through the Bioconda channel, which can be downloaded and installed using the instructions in the Conda documentation.

VOE development

The overview of the VOE development workflow in Fig. 1 summarizes the following process. To begin the process, genomic SARS-CoV-2 SRA accession numbers were selected from the NCBI SRA database filtered with paired library layout, illumina platform and 31/5/2022–31/5/2023, creating input files. Using the NCBI SRA toolkit, 1,961 SRA files were split into paired-end FASTQ files. To create intermediate files, the paired-end FASTQ files were referenced and assembled using the HaVoc pipeline, which includes a read trimming pre-processing tool (Truong Nguyen et al., 2021). Subsequently, 1,200 assembled SRA files were identified as belonging to the SARS-CoV-2 Omicron lineage, using Pangolin version 4.3 (Rambaut et al., 2020). The Omicron lineage accession number was selected. The Omicron SRA files in csv format were split into paired-end FASTQ files using the NCBI SRA toolkit (https://trace.ncbi.nlm.nih.gov/Traces/sra/sra.cgi?view=software). All sequence reads in the SRA files were trimmed to remove low-quality reads and bases, and adapter sequences removed using Fastp (Chen et al., 2018) with a cut-off from the SARS-CoV-2 Freebayes pipeline.

Before the next step, the trimmed sequence reads were validated. The quality of the trimmed reads of Omicron SRA files were assessed using the following four statuses of FastQC version 0.12.1 (Simon, 2010): (1) Per base sequence quality (median value of each base greater than 25), (2) per sequence quality (median quality greater than 27), (3) per base N content (N base less than 5% at each read position) and (4) adapter content (adapter sequences at each position less than 5% of all reads) (Simon, 2010). The output FastQC files were combined with MultiQC v1.4 (Ewels et al., 2016). Files that failed the quality assessment were removed for the variant calling step.

To create the BCFtools variant database, selected fastp (trimmed fastq) files were processed using Bowtie2 Aligner (Langmead & Salzberg, 2012) and BCFtools variant caller (Li, 2011). All reads in fastp files were aligned against Wuhan-Hu-1 (NCBI accession: NC_045512.2) using Bowtie2 version 2.3.4.1 (Langmead & Salzberg, 2012). The SAM files were converted to BAM files and the BAM files were sorted and indexed using SAMtools (Li et al., 2009). The BAM files were accessed by BCFtools version 1.17 (Li, 2011) as variants (read depth > 100) (Lythgoe et al., 2021). All VCF files from the variant calling step were merged with Jacquard version 1.1.5 (Farkas et al., 2021) and annotated with SnpEff version 5.0e (Cingolani et al., 2012). The merged, annotated file was called BCFtools variant database. For the SARS-CoV-2 Freebayes variant database, all selected SRA files were variants called one at a time from the SARS-CoV-2 Freebayes pipeline.

The raw variant databases were then ready to be modified by VOE. The modified BCFtools and SARS-CoV-2 Freebayes variant databases were presented in the format shown in Table 1. Headers and uninteresting information were deleted. POS and #CHROM were swapped. In the INFO column, the variant was annotated using the format of SnpEff version 5.0e. Each database was used as a variant database for analysis of variants on epitopes with VOE. The first database was the BCFtools variant database and the second was the SARS-CoV-2 Freebayes variant database.

VOE was coded with a Python script that locates variants on the epitope and calculates the sensitivity of these variants. The workflow for VOE is shown in Fig. 2. It uses BLAST (Camacho et al., 2009) tools that are publicly available in Bioconda on Unix/Linux platforms.

First, the nucleotide position of the epitope sequence from tBLASTn is generated using the default E-value (E-value <  10) (Camacho et al., 2009) (%identity = 100) and compared with the CDS in the SARS -CoV-2 database. The first hit from the tBLASTn result is selected. The reference database was downloaded from CDS nucleotide Wuhan-Hu-1 (NCBI Accession: NC_045512.2).

Second, all nucleotide positions of the epitope in the SARS-CoV-2 genome are calculated using the database according to the formula: nucleotide position in the gene from the tBLASTn result + nucleotide position in the genome from the database + 1.

Third, all variants on an epitope with HIGH or MODERATE effects (Cingolani et al., 2012) are included in the information of each variant from the variant database because both effects are assumed to change in amino acid sequence. All values in each column of the SRA accession number from the variant database are stored in a dictionary variable, such as {SRR: [0,0,0,1]} (0 means a variant item was not found in the SRA escrow number, 1 means a variant was found).

Last, using the dictionary variables, sensitivity is calculated based on the hypothesis that if a HIGH or MODERATE impact variant is found on the epitope, a false negative result (FN) follows, whereas if no variant or a synonymous variant, or both, are found on the epitope, a true positive result (TP) follows. If any of the values in the columns is 1, then this SRA results in a FN. Conversely, if all of the values in the columns are 0, then this SRA results in a TP. The sensitivity of the epitope is calculated using the following formula (Trevethan, 2017). (1)${\displaystyle Sensitivity=\frac{TP}{TP+FN}.}$

For example, suppose an epitope is located at nucleotide positions 10–30 and in this position interval, no variant is found for SRRX1, one variant is found for SRRX2, one variant is found for SRRX3, three variants are found for SRRX4, and one variant is found for SRRX5. Then, TP = 1 and FN = 4, and sensitivity is 1/(1+4) or 0.2 and %sensitivity is 20%.

Peptide-protein docking validation

Peptides of epitopes were modeled with PEP-FOLD 3 Server (Lamiable et al., 2016). Peptide-protein docking of epitopes with lower 90% sensitivity and the most frequently occurring variant against specific proteins were analyzes using AutoDocktools version 1.5.7 (Morris et al., 2009). An increase in binding energy for mutant docking will make less compact binding of the epitope.

Identification of the SARS-CoV-2 lineage and trimming

The sequence reads of selected 1,791 SRA files were downloaded and assembled as 1,659 consensus files, which were identified with 1,023 Omicron accession numbers. The quality of reads of the accession numbers were assessed, giving 1,011 accepted accession numbers.

VOE variant analysis and sensitivity calculation

Seven epitopes including KLNDLCFTNV, ITLCFTLKRK, GGDGKMKD, RVQPTES, APGQTGK, DSKVGGNYN, and LKPFERD were obtained to evaluate the performance of the proposed VOE analysis tool. The epitopes KLNDLCFTNV, ITLCFTLKRK, RVQPTES, APGQTGK, DSKVGGNYN, and LKPFERD are highly antigenic for the development of an epitope-based peptide vaccine and the epitope GGDGKMKD was identified as ideal for serodiagnostic testing.

The tBLASTn results of all epitopes from the BCFtools variant database appear in Table 2.

The first epitope KLNDLCFTNV is located on the S gene at nucleotide positions 1,156 to 1,185 (amino acid positions 386 to 395 on the S gene) with an E-value of 4.81 × 10⁻⁴. The second epitope ITLCFTLKRK is located at nucleotide positions 328 to 357 on the ORF7a gene (amino acid positions 110 to 119) with an E-value of 0.005. The third epitope RVQPTES is located at nucleotide positions 955 to 975 on the S gene (amino acid positions 319 to 325) with an E-value of 0.025. The fourth epitope APGQTGK is located at nucleotide positions 1,231 to 1,251 on the S gene (amino acid positions 411 to 417) with an E-value of 0.028. The fifth epitope DSKVGGNYN is located at nucleotide positions 1,324 to 1,350 on the S gene (amino acid positions 442 to 450) with an E-value of 0.003. The sixth epitope LKPFERD is located at nucleotide positions 1,381 to 1,401 on the S gene (amino acid positions 461 to 467) with an E-value of 0.008. The last epitope GGDGKMKD is located at nucleotide positions 286 to 309 on the N gene (amino acid positions 96 to 103) an E-value of 0.12. The % identity of all epitopes was 100%.

The variant information for all epitopes from the BCFtools variant database appear in Table 3. For the epitopes KLNDLCFTNV (S386-395) and LKPFERD (S461-467) on the S gene and the epitope ITLCFTLKRK (ORF7a110-119) on the ORF7a gene, no variants with two effects (HIGH or MODERATE) were found, resulting in a sensitivity of 100%. In contrast, for the epitopes RVQPTES (S319-325), APGQTGK (S411-417), DSKVGGNYN (S442-450) on the S gene, and GGDGKMKD (N96-103) on the N gene, variants were found with two effects, and one (p.Thr323Ile), three (p.Lys417Thr, p.Lys417Met, p.Lys417Asn), seven (p.Lys444Thr, p.Lys444Asn, p.Val445Leu, p.Val445Ala, p.Gly446Ser, p.Gly446Asp, p.Asn450Asp), and four (p.Gly96Cys, p.Gly96Val, p.Gly97Cys, p.Gly99Cys) missense variants, respectively, giving respective sensitivities of 99.9011%, 28.4866%, 87.7349%, and 99.6044%.

The tBLASTn results of all epitopes from the SARS-CoV-2-Freebayes variant database appear in Table 4. The first epitope KLNDLCFTNV is located at nucleotide positions 1,156 to 1,185 on the S gene (amino acid positions 386 to 395 on the S gene) with an E-value of 4.81 ×10⁻⁴. The second epitope ITLCFTLKRK is located at nucleotide positions 328 to 357 on the ORF7a gene (amino acid positions 110 to 119) with an E-value of 0.005. The third epitope RVQPTES is located at nucleotide positions 955 to 975 on the S gene (amino acid positions 319 to 325) with an E-value of 0.025. The fourth epitope APGQTGK is located at nucleotide positions 1,231 to 1,251 on the S gene (amino acid positions 411 to 417) with an E-value of 0.028. The fifth epitope DSKVGGNYN is located at nucleotide positions 1,324 to 1,350 on the S gene (amino acid positions 442 to 450) with an E-value of 0.003. The sixth epitope LKPFERD is located at nucleotide positions 1,381 to 1,401 on the S gene (amino acid positions 461 to 467) with an E-value of 0.008. The last epitope GGDGKMKD is located at nucleotide positions 286 to 309 on the N gene (amino acid positions 96 to 103) an E-value of 0.12. The % identity of all epitopes was 100%.

The variant information of all epitopes from the SARS-CoV-2-Freebayes variant database was presented in Table 5. For the epitope KLNDLCFTNV (S386-395), variants with two effects, three missense variants (p.Asn388Tyr, p.Asp389Val, and p.Leu390Pro), and one frameshift variant (p.Thr393fs) were found on the S gene. For the epitope ITLCFTLKRK (ORF7a110-119) on the ORF7a variant, five missense variants (p.Ile110Thr, p.Phe114Val, p.Thr115Ile, p.Leu116Ile, and p.Lys119Arg) were found. For the epitope RVQPTES (S319-325) on the S gene variant, one stop-gain (p.Arg319*), one frameshift (p.Pro322fs), and two missense variants (p.Val320Phe, and p.Thr323Ile) were found. For the epitopes APGQTGK (S411-417), DSKVGGNYN (S442:450), and LKPFERD (S461-467) on the S gene variant, seven (p.Gln414Lys, p.Lys417Asn, p.Lys417His, p.Lys417Asn, p.Lys417Thr, p.Lys417Ile, and p.Lys417Asn), 10 (p.Lys444Thr, p.Lys444Asn, p.ValGly445ProSer, p.Gly446Ser, p.Val445Ala, p.ValGly445AlaSer, p.Gly446Ser, p.Gly446Asp, p.Asn448Ile, and p.Asn450Asp) and two missense variants (p.Pro463Leu, and p.Asp467Val) were found, respectively. For the epitope GGDGKMKD (N96-103) on the N gene variant, four missense variants (p.Gly96Cys, p.Gly96Val, p.Gly97Cys, and p.Gly99Cys) were found. Lastly, the respective % sensitivities of all epitopes were 99.5054%, 99.5054%, 99.6044%, 6.8249%, 71.1177%, 99.8022%, and 99.6044%.

Peptide-protein docking validation

Epitopes with low sensitivity (<90%) could be considered moderate to low binders and were further investigated. The 3D structures of the B cell receptor (BCR) (ID: 5IFH) (Rantam et al., 2021) were obtained from the RCSB Protein Data Bank. The wild-type epitope (APGQTGK) and the mutant epitope (APGQTGN) were docked with the BCR shown in Figs. 3A and 3B, respectively. The estimated binding energy was −2.56 kcal/mol for the wild-type docking and −2.02 kcal/mol for the mutant docking. The wild-type epitope (DSKVGGNYN) and mutant epitope (DSKVSGNYN) were docked with the BCR shown in Figs. 4A and 4B, respectively. The estimated binding energy was −1.69 kcal/mol for wild-type docking and +3.58 kcal/mol for mutant docking.