{'fraction_non_alphanumeric': 0.029358687747461417, 'fraction_numerical': 0.02663829558003609, 'mean_word_length': 4.772388059701493, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Bioinformatics pipeline unveils genetic variability to synthetic vaccine design for Indian SARS-CoV-2 genomes

In the worrisome scenarios of various waves of SARS-CoV-2 pandemic, a comprehensive bioinformatics pipeline is essential to analyse the virus genomes in order to understand its evolution, thereby identifying mutations as signature SNPs, conserved regions and subsequently to design epitope based synthetic vaccine. We have thus performed multiple sequence alignment of 4996 Indian SARS-CoV-2 genomes as a case study using MAFFT followed by phylogenetic analysis using Nextstrain to identify virus clades. Furthermore, based on the entropy of each genomic coordinate of the aligned sequences, conserved regions are identified. After refinement of the conserved regions, based on its length, one conserved region is identified for which the primers and probes are reported for virus detection. The refined conserved regions are also used to identify T-cell and B-cell epitopes along with their immunogenic and antigenic scores. Such scores are used for selecting the most immunogenic and antigenic epitopes. By executing this pipeline, 40 unique signature SNPs are identified resulting in 23 non-synonymous signature SNPs which provide 28 amino acid changes in protein. On the other hand, 12 conserved regions are selected based on refinement criteria out of which one is selected as the potential target for virus detection. Additionally, 22 MHC-I and 21 MHC-II restricted T-cell epitopes with 10 unique HLA alleles each and 17 B-cell epitopes are obtained for 12 conserved regions. All the results are validated both quantitatively and qualitatively which show that from genetic variability to synthetic vaccine design, the proposed pipeline can be used effectively to combat SARS-CoV-2.


Introduction
More than two years ago, SARS-CoV-2 put a massive halt to the freedom of human movement due to its high transmission rates [1]. Early study established the fact that SARS-CoV-2 virus is highly similar to that of the SARS-CoV-1 (95%-100%) [2]. In April 2021, India registered its second sudden surge in official cases with the Delta (B.1.617.2) variant. In late 2021, the third wave hit the country which was led by Omicron. Though, India is pushing towards a very large vaccination drive, concerns over the efficacy of the vaccine for such aggressive mutations are also increasing. Meanwhile, India is not the only country which has witnessed the new mutation strain of the evolving virus, variants in South Africa (501Y.V2) [3], United Kingdom (B.1.1.7) [4], Japan (E484K) [5], Brazil (P.1) [5] are also making their rounds. The latest variant to join the bandwagon is Omicron (B.1.1.529). Although, previously it was suggested that such mutants are not going to affect the effectiveness of the vaccines currently in use, the emergence of Omicron has changed the equation. Moreover, new variants can affect the diagnosing procedure such as primer identification or antibody binding in RT-PCR. Ascoli [6] also suggested that a mutation in the Spike region of SARS-CoV-2 may affect the diagnosing procedure with greatest impact along with the increasing infection rates, transmissibility or even impacting people of younger age.
In the current scenario, it is an important and urgent task to study the frequently occurring mutations within the virus. In this regard, Yuan et al. [7] have analysed 11,183 SARS-CoV-2 genome from around the globe to identify the SNPs and critical SNPs with specific high mutation frequency along with the geographical pattern analysis. Further, they have found 74 non-synonymous and 43 synonymous mutations. Most importantly they have identified Nucleocapsid (N) as the gene with the highest mutational frequency changes. This directly undermines the claim of Ascoli [6] that Nucleocapsid can be targeted for the diagnosing purposes as N gene undergoes very less mutations or is mostly conserved. Hence, it is important to take a closer look how SARS-CoV-2 is evolving over time. Moreover, Tang et al. [8] have found new developing variations on the receptor binding sites of Spike gene of SARS-CoV-2 in the form of S and L lineages. Here, S and L lineages are defined by two tightly linked SNPs at positions 8,782 (orf1ab:T8517C, synonymous) and 28,144 (ORF8: C251T, S84L) which might affect the virus pathogenesis. Phylogenetic analysis done by Maitra et al. [9] revealed the signature mutations such as C14408T in RdRp along with A23403G change in Spike protein majorly forming A2a clade within 9 Indian sequences. Further, they have also reported a triplet based mutation in N gene 2881-3 GGG/AAC which might affect the miRNAs bindings to original sequences. Genome analysis by Saha et al. [10] for 72 different countries has shown multiple unique mutation points in the form of substitution, deletion, insertion and SNPs in each country, resulting in 7209, 11700, 119 and 53 mutations respectively. Further, they have identified 11 SNPs which are unique to India, the most frequent being T1198K, A97V, T315N and P13L mutation points in NSP3, RdRp, Spike and ORF8. Therefore, it has become more important than ever to constantly monitor the continuous evolving virus in order to take up proper measures to battle the contagious virus. Study conducted by Nagy et al. [11] identified genomic alterations and the association of each mutation and outcome. As a result, they have found 3733 mutation points related to mild outcome in ORF8, NSP6, ORF3a, NSP4 and Nucleocapsid genes whereas the mutations in Spike glycoprotein, RNA polymearse, ORF3a, NSP3, ORF6 and N provided inferior outcome. Also, severe outcomes are associated to the mutations in ORF3a and NSP7 proteins. Thus, mutations are important in the significant genes such as Spike, N etc. and such mutations may even lead to a false diagnosis in RT-PCR testing. Hence, it is also important to extract the conserved regions in a genomic sequence for more effective diagnosis. In this regard, [10] have identified a conserved region in NSP6 gene as a potential target for SARS-CoV-2 detection using RT-PCR.
On the other hand, alteration in the RNA virus can lead to vaccine failures as was noticed in the case of Influenza virus in 2013-14 [12]. Hence, to fight against a highly evolving virus like SARS-CoV-2, it is important to have stable vaccine. In this regard, Ghosh et al. [13] have performed a genome-wide analysis of 10644 SARS-CoV-2 sequences to identify the conserved regions in a virus genome, followed by which they have proposed epitope based vaccine design targeting the T-cell and B-cell epitopes. Another study conducted by Ghosh et al. [14] for identifying the conserved regions specifically focussed on 566-Indian SARS-CoV-2 sequences by considering four different multiple sequence alignment techniques. In both the studies most immunogenic and antigenic epitopes were derived from various coded proteins of the virus which can be targeted for synthetic vaccine design. Alam et al. [15] targeted the Spike glycoprotein to propose non-allergic, highly antigenic and non-mutant synthetic vaccine design targeting Thymus cell (T-cell) and bone marrow. Rahman et al. [16] targeted 3 important genes viz Spike, Membrane and Envelope for multi-epitope-based vaccine design for SARS-CoV-2 with a 90% population coverage. Also, immune simulation suggested a significant increase in primary immune response with increased IgM and secondary immune response with increased IgG1 and IgG2 along with increased proliferation of T-helper cells with increased cytokines. Another study [17] targeted heptad repeats 1 and 2 (HR1 and HR2) in the Spike protein for peptide design using molecular dynamics simulation between the fusion of the viral membrane with the host cell membrane. This eventually limited the spread of the virus in the host cells. Vashi et al. [18] predicted 24 potential epitope fragments of which 20 were on the surface of Spike protein (S protein) and were considered to be helpful for designing potential immunogenic peptide based vaccines.
Motivated by the literature and looking at the sudden surge of SARS-CoV-2 in India, a comprehensive bioinformatics pipeline is proposed in this work to analyse the virus genomes for understanding its evolution for identifying mutations as signature SNPs, conserved regions and subsequently to design epitope based synthetic vaccine. In this regard, we have performed multiple sequence alignment of 4996 Indian SARS-CoV-2 sequences as a case study using MAFFT followed by phylogenetic analysis of the aligned sequences using Nextstrain. As a result, the sequences are found to be distributed in 5 clades, viz 19A, 19B, 20A, 20B and 20C. Thereafter, from the aligned sequences, mutation points as SNPs are identified in each clade. Subsequently, top 10 signature SNPs based on their frequency are identified in each clade resulting in a total of 50 such SNPs. Out of 50 signature SNPs, 40 unique signature SNPs are identified resulting in 23 non-synonymous signature SNPs which gives 28 amino acid changes in protein which are visualised in protein structures as well. Furthermore, the sequence and structural homologybased prediction along with the protein structural stability of the amino acid changes for such SNPs are evaluated using PROVEAN, PolyPhen 2.0 and I-Mutant 2.0 in order to judge the characteristics of the identified clades. As a consequence, A97V in RdRp in 19A, V354L in Nucleocapsid in 19B, Q57H in Nucleocapsid in 20A, R203M in Nucleocapsid in 20B while T85I in NSP2 and Q57H in ORF3a in 20C are the unique amino acid changes which are responsible for defining each clade as they are all deleterious and unstable and are also responsible for decreasing the protein structural stability. Moreover, based on the entropy of each genomic coordinate of the aligned sequences, conserved regions are identified. Conserved regions are such places in genomic sequences for which the corresponding protein sequences remain unchanged. These conserved regions are then filtered based on the criteria that their lengths are greater than or equal to 125nt and their BLAST specificity score is equal to 100% resulting in 12 conserved regions belonging to NSP2, NSP8, NSP10, RdRp, Exon, Spike glycoprotein, ORF3a and ORF7a proteins. Based on its length, one conserved region as potential target is identified in the NSP10 gene for which the primers and probes are reported as well. Such primers and probes can be used for detecting SARS-CoV-2 virus. The 12 conserved regions are also used to identify the T-cell and B-cell epitopes along with their immunogenic and antigenic scores. Using such scores, most immunogenic and antigenic epitopes are selected for the 12 conserved regions thereby identifying 23 MHC-I and 22 MHC-II restricted T-cell epitopes with 10 unique HLA alleles each and 17 B-cell epitopes. Finally, the binding conformation of the MHC-I and MHC-II restricted T-cell epitopes with respect to HLA alleles are shown to judge their relevance. Also, the physico-chemical properties of the epitopes are reported along with structural properties using Ramchandran plots, ERRAT score and Z-Scores. Thus, based on the comprehensive bioinformatics pipeline, the main contributions of this work can be summarised as: (a) phylogenetic analysis in Nextstrain to identify virus clades, (b) identification of SNPs in the aligned sequences, (c) based on frequency, top 10 signature SNPs identification in each virus clade, (d) identification of conserved regions and based on length selecting one such region as potential target for reporting the corresponding primers and probes to detect SARS-CoV-2 and (e) identification of T-cell and Bcell epitopes for peptide based synthetic vaccine design.

Material and Methods
In this section, the details of data collection and the preparation are elucidated which is followed by a brief discussion on the pipeline of the workflow that has been considered in this work.

Data Collection and Preparation
The reference sequence of SARS-CoV-2 virus (NC_045512.2) is collected from National Center for Biotechnology Information (NCBI) 2 while 4996 complete or near complete Indian SARS-CoV-2 genomes are collected from Global Initiative on Sharing All Influenza Data (GISAID) 3 in fasta format. The 4996 SARS-CoV-2 sequences are mostly distributed from January 2020 to January 2021. These sequences are then aligned to find the conserved regions. The coded protein corresponding to each conserved region is extracted as well. Further, to map the protein sequences and changes in the amino acid, protein PDB are collected from Zhang Lab 4 which are then used to model and identify the structural changes. All these analyses are executed on High Performance Computing (HPC) facility of NITTTR, Kolkata while the amino acid changes are checked in MATLAB R2019b.

Pipeline of the work
The pipeline of this work is provided in Fig. 1. In this work, a comprehensive bioinformatics pipeline is proposed which encompasses identifying mutation points as SNPs, conserved regions and finally design of epitope based synthetic vaccine. To achieve these goals, in the first phase of the pipeline, multiple sequence alignment of 4996 Indian SARS-CoV-2 genomes as a case study using MAFFT [19] is carried out followed by the phylogenetic analyses using Nextstrain [20]. As MAFFT uses fast fourier transform, it outperforms all the other alignment techniques. On the other hand, analysis of the evolution and spread of pathogens is done using Nextstrain by considering phylogenomic and phylogeographic data. The spread and evolution of virus genomes can be visualised at nextstrain.org using auspice. By using this tool, the evolution and geographic distribution of SARS-CoV-2 genomes are visualised by creating the metadata in our High Performance Computing environment. Once the identification of the virus clades are performed using Nextstrain, clade specific aligned sequences are used to identify mutation points as substitutions especially SNPs in each clade. Henceforth, codon table is used to identify the amino acid changes in the virus proteins corresponding to the SNPs. Thereafter, based on their frequency in the virus genome, top 10 signature SNPs are identified in each clade. Please note that the amino acid changes in the SNPs can be either synonymous or non-synonymous. Furthermore, amino acid changes in the non-synonymous SNPs are visualised in the protein structures and they are used to evaluate their functional characteristics as well.
The second phase of the pipeline entails identification of conserved Regions (CnRs) in the aligned sequences using entropy (N ) which can be computed as: where F y x represents the frequency of each residue x occurring at position y and 5 represents the four possible residues as nucleotides plus gap. To identify the conserved regions, a minimum segment length of 15 is considered with maximum average entropy as 0.2 along with a maximum entropy per position of 0.2 as well without any gaps. All these values are taken after following the literature. Thereafter, refinement criteria for the conserved regions are adopted based on the criteria that their lengths are ⩾125nt and their BLAST specificity score as query coverage is equal to 100%. Subsequently, based on its length, a particular conserved region is considered as potential target which is then used to identify primers and probes using Primer-BLAST 5 for SARS-CoV-2 detection.
In the final phase of the pipeline, T-cell and B-cell epitopes along with their immunogenic and antigenic scores are predicted for the refined CnRs using IEDB 6 and ABCPred 7 respectively. For such MHC-I and MHC-II restricted T-cell epitopes, predictions are carried out using IEDB recommended NetMHCPan EL 4.1 8 and Consensus Approach 9 [21] respectively while ABCPred [22] is used for B-cell epitope prediction. Thereafter, by using these predicted epitopes, antigenic scores are evaluated by VaxiJen 2.0 10 while the validation of the identified T-cell epitopes is carried out by studying their conformational 2D non-covalent structures using LigPlot+ [23]. For the verification of the predicted Bcell epitopes, BepiPred 2.0 11 [24] server is used. Allergen and toxicity properties of the epitopes are evaluated using AllerTop 2.0 12 and Tox-inPred 13 respectively. The physico-chemical properties are also evaluated using ToxinPred. Moreover, docking of all the T-cell epitopes are performed using AutoDock Vina [25] and their structural properties are reported using Ramachandran Plot [26], ERRAT score [27] and Verify 3D [28] using SAVES 6.0 14 . Finally, Z-Score evaluation is performed using ProSA [29].

Phylogenetic analysis and Signature SNPs in each clade
To achieve the first step of the bioinformatics pipeline, multiple sequence alignment of 4996 Indian SARS-CoV-2 genomes is performed using MAFFT followed by phylogenetic analysis with the help of Nextstrain. This phylogenetic analysis results in 5 clades viz. 19A, 19B, 20A, 20B and 20C. Thereafter, mutation points as substitutions specifically SNPs are identified in each clade resulting in 708, 161, 3308, 3235 and 47 SNPs for 479, 88, 2486, 1925 and 18 sequences respectively in 19A,  Table S1. The resultant phylogenetic trees in radial and rectangular views are shown in Fig. 2(a) and (c) while the clade wise geographical distribution of the 4996 sequences is shown in Fig. 2 The clade wise evolution of the sequences for each month of each Indian state is shown in the form of pie charts in supplementary Table S2 while the month wise evolution of such sequences for each clade is reported in supplementary Table S3. The corresponding colour representation for the five major clades and the months are provided in supplementary Figure S1. Moreover, the entropy values for the nucleotide changes and coding regions of the SARS-CoV-2 genome are shown respectively in Fig. 2(d) and (e). It is to be noted that for some sequences, the state name is not mentioned in the GISAID database. Thus, they are aggregated under the state name 'India'.
Once the SNPs are determined for each clade, top 10 SNPs based on their frequency viz. signature SNPs are identified in each clade, thereby resulting in 50 signature SNPs as reported in Table 1 and visualised in Fig. 2(f). In unsupervised learning, feature selection is a very crucial task. In this work, frequency of a SNP is considered to be the feature selection criterion. For example, G11083A and G11083T with a frequency of 425 is the top signature SNP in clade 19A while for 19B, T28144C having frequency of 87 is the top signature SNP. Subsequently, 40 unique SNPs are identified which results in 23 non-synonymous signature SNPs with 28 corresponding amino acid changes. The common signature SNPs in the five clades are visualised using Venn diagram in Fig. 2(g). It is evident from the figure that the clades do not have any common SNPs, thereby confirming the fact that signature SNPs are indeed the defining features of a clade. Moreover, the amino acid changes are visualised in Fig. 3 as well. Please note that 27 amino acid changes are visualised in Fig. 3 as opposed to 28 reported changes; the discarded change is E110* in ORF8 as this amino acid change leads to a stop codon. Also, sequence and structure-based homology prediction of the amino acid changes for the non-synonymous SNPs are reported in Table 2, the details of which are discussed in Discussion section. All the  Table S1.

Selection of CnRs
For the next phase of this study, we have obtained 473 conserved regions (CnRs) which are then mapped to the 11 coding regions of SARS-CoV-2; ORF1ab, Spike, ORF3a, Envelope, Membrane, ORF6, ORF7a, ORF7b, ORF8, Nucleocapsid and ORF10. For each CnR, the corresponding protein sequence is taken according to the reading frame it is associated with. For example, protein sequence of CnR in Spike region is taken from Frame 2 while that belonging to Envelope and Membrane are taken from Frames 1 and 3 respectively. These 473 conserved regions are then filtered based on the criteria that the length of the CnR should be greater than or equal to 125nt and the their BLAST specificity score as query coverage is equal to 100%. As a result, we have obtained 12 such regions as reported in Table 3. The table also shows the corresponding protein sequences for the conserved regions along with their length, BLAST specificity score, percent of BLAST specificity score as query coverage, coding regions, starting and ending coordinates, length of coding regions and the coded proteins. These CnRs belong to coding regions which code NSP2, NSP8, NSP10, RdRp, Exon, Spike glycoprotein, ORF3a protein and ORF7a protein. The details of all the initial and filtered CnRs are provided in the supplementary as an excel file. Also, based on its length, among these CnRs, one CnR is then chosen as the target for the detection of SARS-CoV-2. Moreover, the protein sequences of these CnRs are used to identify the MHC-I and MHC-II restricted T-cell and B-cell epitopes.

Identification of Conserved Region as Target and associated Primers and Probes
Among the 12 CnRs identified, the CnR with the largest length of 247nt is considered to be a potential target. This CnR belongs to ORF1ab region, specifically NSP10 gene which is shown in Table 4. With a Nucleotide BLAST score of 457 and BLAST specificity score as query coverage is equal to 100%, the global stability of this CnR as a global target is confirmed. The structure of the NSP10 gene as shown in Table 4 is taken from ZhangLab in the form of a PDB file and the CnR as target is highlighted in red. Using this conserved region, 10 primers and probes are identified from Primer-BLAST and reported in Table 5 and shown in Fig. 2(h). The table reports both the forward and the reverse primers. Moreover, high GC scores (45%-53%) of the identified primers suggest that the identified primers and probes can be used in RT-PCR for SARS-CoV-2 detection in order to correctly diagnose COVID-19 patients. Therefore, the target region of NSP10 gene can be considered as a confirmatory assay. It is to be noted that based on its adhesive properties, Ong et al. [30] have predicted NSP10 as a possible vaccine candidate.

Identification of T-cell Epitopes
To achieve the final phase of the pipeline, design of epitope based synthetic vaccine is carried out. To predict the epitopes from the 12 CnRs, the corresponding protein sequences are fed to the various tools as inputs. For the prediction of MHC-I restricted T-cell epitopes, IEDB recommended NetMHCPan EL 4.1 [31] is considered targeting 27 unique HLA alleles. For each CnR, this resulted in the selection of 5 best HLA allele binder epitopes based on their immunogenic scores. Thereafter, these best binders are provided as input to VaxiJen 2.0 [32] server for antigenic score prediction [31] with a cut-off score of 0.4. Any epitope beyond this cut-off are considered to be antigenic. Therefore, a total of 60 epitopes, each of length 9-10 mer, are obtained along with their immunogenic and antigenic scores. From each of the 12 CnRs, the most immunogenic and antigenic MHC-I restricted T-cell epitopes are identified resulting in 22 such epitopes and reported in Table 6. With a score of 0.99, the most immunogenic epitopes are SEVGPEHSL, DTDFVNEFY and QEYADVFHLY bounded to HLA-B*40:01, HLA-A*01:01 and HLA-B*44:03 alleles respectively belonging to NSP2 and RdRp coded proteins. On the other hand, with a score of 1.43, HPNPKGFCDL is the most antigenic epitope belonging to NSP10 coded Table 2 Sequence and structural homology-based prediction for non-synonymous signature SNPs along with their protein structural stability.   Similarly, MHC-II restricted T-cell epitopes are predicted using IEDB recommended consensus approach targeting a different set of 27 unique HLA alleles resulting in 60 epitopes, each of length 15 mer. Subsequently, the most immunogenic and antigenic MHC-II restricted T-cell epitopes are identified for the 12 CnRs which resulted in 21 such epitopes as reported in Table 6. It is to be noted that a MHC-II restricted Tcell epitope with a low immunogenic score is a better vaccine candidate. Thus, with a score of 0.02, NEFYAYLRKHFSMMI belonging to RdRp coded protein and bounded to HLA-DRB1*11:01 allele is the most immunogenic epitope while the most antigenic epitope is GCVPLNIIPLTTAAK belonging to NSP8 coded protein and bounded to HLA-DRB1*08:02 allele. All the 60 MHC-I and MHC-II restricted T-cell epitopes along with their HLA alleles are provided in the supplementary as an excel file and the corresponding link is provided in Table S1.

Identification of B-cell Epitopes
Epitope designing consists of both T-cell as well as B-cell epitopes; the latter one is particularly important for antigen production against a virus. In this regard, ABCPred is used for the prediction of B-cell epitopes where a threshold of 0.5 is maintained to consider the epitopes beyond this threshold to be immunogenic. With a cut-off value of 0.4, VaxiJen 2.0 server is used to evaluate the antigenic scores of the epitopes. Thus, we have identified 50 linear B-cell epitopes, each of length 16 mer, for the 12 CnRs, among which 17 are selected to be the most immunogenic and antigenic as shown in Table 6. These epitopes are also verified with the help of BepiPred 2.0 server and their corresponding graphical analysis is shown in supplementary Figure S2 where the red line represents the threshold which is set to 0.35 and the total green and yellow regions indicate a protein sequence. The most immunogenic and antigenic B-cell epitopes as reported in Table 6 are respectively VVKIYC-PACHNSEVGP belonging to NSP2 coded protein and PNPKGFCDLKGKYVQI belonging to NSP10 coded protein. Their corresponding graphical representations are provided in supplementary Figure S2 (a) and (c) respectively. All the 50 B-cell epitopes are provided the supplementary as an excel file and the corresponding link is provided in Table S1.
Additionally, in Table 7 we have provided a summarised list of all the epitopes belonging to these 12 CnRs along with their allergic and toxicity characteristics predicted using AllerTOP 2.0 15 and ToxinPred 16 where 12, 6 and 8 allergic MHC-I, MHC-II T-cell and B-cell epitopes are identified respectively while only 1 and 5 epitopes in MHC-I restricted Tcell and B-cell epitopes are found to be toxic. The 3D structures of the epitopes summarised in Table 7 are further highlighted in Fig. 4 using ChimeraX. For better understandability, the identified epitopes are underlined in supplementary Figure S3.

Discussion
Since its emergence in Wuhan, China, SARS-CoV-2 has spread very rapidly around the world resulting in a global pandemic. Though the vaccination process has started, the number of COVID affected patients is still quite large. The waves of COVID-19 pandemic are a huge threat to the human population. In this regard, it is important to develop a bioinformatics pipeline in order to conduct in-depth analysis of SARS-CoV-2 genomes in every one or two months for next four to five years to know the evolution, genetic variability, virus strains and conserved regions, thereby to use such information for proper vaccine. Moreover, the mutated variants found in India are also a major concern of the researchers. Thus, identification of virus strains is very essential in today's scenario. Moreover, vaccine is the only ray of hope in this dire situation, Table 4 Targeted Conserved Region in SARS-CoV-2 Genome and its corresponding protein sequence in NSP10 which is highlighted by red colour in NSP10 gene. Indian SARS-CoV-2 genomes which has resulted in the identification of five clades and subsequently 10 signature SNPs in each clade. Also, based on entropy, conserved regions are identified for the aligned sequences and primers and probes are identified as well for SARS-CoV-2 detection. Furthermore, we have identified T-cell and B-cell epitopes for the development of vaccines. Structural changes in amino acid residues can often result in changes in the protein translations which is conducive to functional instability of the proteins. In this regard, sequence and structural homology-based prediction of the amino acid changes in the non-synonymous signature SNPs along with their protein stability for the 4996 sequences are reported in Table 2 using PROVEAN (Protein Variation Effect Analyser) [33], PolyPhen-2 (Polymorphism Phenotyping) [34] and I-Mutant 2.0 [35] to judge the characteristics of the identified clades. PROVEAN 17 works with sequence based prediction algorithm while Polyphen-2 18 uses prediction based on sequence, structural and phylogenetic information of a SNP. I-Mutant 2.0 19 uses support vector machine (SVM) for  Table 6 List of most Immunogenic and Antigenic Epitopes for MHC-I, MHC-II restricted T-cell and B-cell Epitopes for 12 CnRs. *I.S.-Immunogenic Score; A.S.-Antigenic Score. the automatic prediction of protein stability changes for SNPs. PRO-VEAN and PolyPhen-2 are used to find the deleterious or damaging nonsynonymous SNPs. The threshold value of PROVEAN is set to − 2.5. If the PROVEAN score of a SNP is le this threshold, the corresponding nonsynonymous mutation is deleterious. For Polyphen-2, this range is between 0 to 1. If the score is closer to 1, mutations are more confidently considered to be damaging. As reported in Table 2 by considering the consensus of PROVEAN and Polyphen-2, out of the 28 unique amino acid changes, 8 unique changes are deleterious and damaging. Moreover, protein stability is important for considering the functional and structural activity of a protein. Any change in protein stability may cause degradation of proteins. The protein stabilities for the nonsynonymous signature SNPs are determined using I-Mutant 2.0. The changes in the protein stability in I-Mutant 2.0 tool is predicted using free energy change values (DDG). A decrease in protein stability is indicated by a zero or a negative value of DDG. Table 2 shows that out of the 8 unique changes, 5 unique changes show a decrease in the stability of the protein structures. As a consequence, A97V in RdRp in 19A, V354L in Nucleocapsid in 19B, Q57H in Nucleocapsid in 20A, R203M in Nucleocapsid in 20B while T85I in NSP2 and Q57H in ORF3a in 20C are the unique amino acid changes which are responsible for defining each clade as they are all deleterious and unstable as well as decrease the stability of the protein structures. All of them are marked in bold in Table 2.
Physico-chemical properties are considered to show the significance of the epitopes as reported in this paper. For each property, the physicochemical values lie between 0 and 1. The physico-chemical properties for MHC-I, MHC-II restricted T-cells and B-cell epitopes belonging to the 12 CnRs are reported in Supplementary Tables S4, S5 and S6 respectively. As reported in Table S4,  For further validation, the conformational 2D non-covalent structures of the MHC-I and MHC-II restricted T-cell epitopes are studied using LigPlot+. Also, it is very important to study their structural characteristics such as binding conformation. Hence, to identify the stable binding interactions, molecular docking of the MHC-I and MHC-II restricted T-cell epitopes are evaluated using AutoDock Vina. For the same, first the 3D structures of the epitopes are prepared with the help of the build structure function of Chimera 1.14 along with the retrieval of the crystal structures of the HLA alleles in PDB format from RCSB Protein Data Bank. To identify the binding energy at the binding groove in the HLA allele, grid search space is set to (60,60,60) with centre of grid at (0,0,0) for X,Y Z coordinate with a spacing parameter of 0.964. The best is selected by higher binding affinity i.e. lowest docking score generated through Autodock Vina. Also, we have used DOE-MBI services such as PROCHECK, ERRAT, Verify3D for Ramchandran Plot, structure quality and 3D structure verification respectively. The results of the docking analysis along with Z-score, respective PDB ID 20 , total energy of the 3D complex, van der Waals energy and electric energy of each complex are reported in Tables 8,9 Table 7 Summary of the most Immunogenic and Antigenic Epitopes along with the Allergic and Toxicity values.  It is to be noted that in our previous works [13,14], with a refinement criteria of 60nt, respectively 17 and 23 conserved regions were identified with 30, 24 and 21 and 34, 37 and 29 best immunogenic and antigenic MHC-I and MHC-II T-cell and B-cell epitopes. These experiments were conducted for SARS-CoV-2 sequences till July 2020. As the virus is constantly evolving, a more recent analysis is needed to understand the evolution of the epitopes. Therefore, this work which uses sequences till January 2021 is very relevant in current scenario of constant virus mutation.

Conclusion
In the past two years, India has witnessed different surges of COVID-19 cases. Hence, it is important to provide a comprehensive bioinformatics pipeline to understand the virus evolution for identifying the mutation points as SNPs, conserved regions and design potential candidates for vaccine design. In this regard, initially, multiple sequence  Table 9 Docking and Z-scores of most Immunogenic and Antigenic MHC-II restricted T-cell epitopes for 12 CnRs.     alignment of 4996 Indian SARS-CoV-2 genomes as a case study are carried out using MAFFT followed by phylogenetic analysis by Nextstrain to identify virus clades, resulting in 5 virus clades; 19A, 19B, 20A, 20B and 20C. Thereafter, mutation points as SNPs are identified in each clade from which top 10 signature SNPs are further identified based on their frequency in each clade. 40 unique signature SNPs are thus identified from the total 50 signature SNPs resulting in 23 non-synonymous signature SNPs which provides 28 amino acid changes in protein. These changes are visualised in their respective protein structure as well. The sequence and structural homology-based prediction of the nonsynonymous signature SNPs along with their protein structural stability are evaluated to judge the characteristics of the identified clades. 40 unique signature SNPs are thus identified from the total 50 signature SNPs resulting in 23 non-synonymous signature SNPs which provide 28 amino acid changes in protein. As a consequence, A97V in RdRp in 19A, V354L in Nucleocapsid in 19B, Q57H in Nucleocapsid in 20A, R203M in Nucleocapsid in 20B while T85I in NSP2 and Q57H in ORF3a in 20C are the unique amino acid changes which are responsible for defining each clade as they are all deleterious and unstable as well as they decrease the protein structural stability. Furthermore, based on the entropy of each genomic coordinate of the aligned sequences, 473 conserved regions are identified which are then refined based on the criteria that their lengths are greater than 125nt and their BLAST specificity score as query coverage is equal to 100%. This refinement results in 12 conserved regions belonging to NSP2, NSP8, NSP10, RdRp, Exon, Spike glycoprotein, ORF3a and ORF7a proteins. Based on length, one conserved region belonging to NSP10 gene is considered to be the potential target for which the corresponding primers and probes are reported for SARS-CoV-2 detection. The 12 conserved regions are then used to identify the T-cell and B-cell epitopes along with their immunogenic and antigenic scores. Such scores are then used to select the most immunogenic and antigenic T-cell and B-cell epitopes resulting in 22 MHC-I and 21 MHC-II restricted T-cell epitopes with 10 unique HLA alleles each and 17 B-cell epitopes. Finally, the relevance of these epitopes are validated by showing the binding conformation of the MHC-I and MHC-II restricted T-cell epitopes with respect to HLA alleles. Also, the physico-chemical properties of the epitopes are reported along with the structural properties using Ramchandran plot, ERRAT scores and Z-Scores. Hence, from genetic variability to synthetic pipeline, a comprehensive bioinformatics pipeline is presented in this study to fight against SARS-CoV-2.

Ethics approval and consent to participate
The ethical approval or individual consent was not applicable.

Availability of data and materials
The aligned 4996 Indian SARS-CoV-2 genomes with the reference sequence and the final results of this work are available at 'http://www. nitttrkol.ac.in/indrajit/projects/COVID-Pipeline-5K/". Moreover, the SARS-CoV-2 genomes used in this work are publicly available at GISAID database..

Consent for publication
Not applicable.

Funding
This work has been partially supported by CRG short term research grant on COVID-19 (CVD/2020/000991) from Science and Engineering Research Board (SERB), Department of Science and Technology, Govt. of India.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.