Jiaqi Wu1 and So Nakagawa2
- The Graduate School of Integrated Sciences for Life, Hiroshima University
- Tokai University School of Medicine
SGV-caller is an open source software under MIT license that can generate the local database of genomic variations of SARS-CoV-2 genomes at amino acid, codon and nucleotide levels. The SGV-caller was designed to handle SARS-CoV-2 genomes and metadata files in the GISAID database. In particular, the update of the GISAID database can be easily tracked. The output files generated by the SGV-Caller are tab-delimited text formats that are easy to read, manage and use. A data quality check is performed on each genome and summary statistics are reported. The data will be stored in a variety of formats that can be easily accessed by R software, bash commands and other programming languages. Based on the outputs, it will be easy to manage the variant information of SARS-CoV-2 genomes and extract the data that meet the user's requirements. SGV-caller supports the use of multiple threads, which could greatly improve the speed of computation. The detail of this software can be also find at: https://doi.org/10.1016/j.heliyon.2025.e42613.
The SGV-Caller software consists of scripts written in the Perl language. Please ensure that you have Perl installed on your computer. In addition, the software requires the multiple alignment software MAFFT (1). It is recommended that you install the latest version of MAFFT. For Ubuntu users, the easiest way to install MAFFT is as follows:
sudo apt-get -y install mafft
If you are using a Mac, you can use brew to install it as follows:
brew install mafft
If you want to install using the source code of the MAFFT suite, you can install it by following the guidance of MAFFT. If this is the case, please set the path to the binary file of MAFFT programs by:
export PATH=$PATH:/path/to/the/directory/where/MAFFT/installed
SGV-Caller also requires xz to be installed in the system. In Linux systems, xz is installed by default. If you are using Mac OS, please install it:
brew install xz
Please download the SGV-caller from https://github.com/wujiaqi06/SGV-caller. You can also run it:
git clone https://github.com/wujiaqi06/SGV-caller.git
The SGV-Caller is a Perl software that does not need to be compiled. First you have to set the path to the SGV-caller directory as follows:
chmod 777 /path/to/the/directory/where/sgv-caller/installed/sgv-caller.pl
To use the SGV-caller, copy the SGV-caller configuration file "sgv-caller.conf" into the directory containing the data to be analysed. Then set the appropriate variables in "sgv-caller.conf".
A "benchmark.pl" is prepared in the "SGV-caller/examples/benchmark/" folder to test the speed of the software (written by Dr Kirill Kryukov in National Institute of Genetics, Japan). You can run it on a Linux system computer by:
perl benchmark.pl
The main purpose of SGV-caller is to generate the local database of genomic variations of SARS-CoV-2 using data downloaded from GISAID database. Variations at nucleotide, amino acid and codon level will be reported. It can also be used for user-defined virus genomes.
First, set the path to the SGV-caller folder and configure the local environment:
export PATH=$PATH:/path/to/the/directory/where/sgv-caller/installed
export LC_ALL=en_US.UTF-8
export LC_CTYPE=en_US.UTF-8
Then copy the configuration file "sgv-caller.conf" to the directory where you hope to generate the database, then modify the values of "sgv-caller.conf" file to control software which calculation you want to perform and its input files, then run this in command line:
cd where/sgv-caller.conf/is/stored
sgv-caller.pl
By default, SGV-caller searches for the configuration file "sgv-caller.conf" in the database directory. If your configuration file has a different name, please run this command instead, then you can use your own configuration file:
sgv-caller.pl -i your_own_configuration_file
Version check for SGV-caller can be done by:
sgv-caller.pl -v
Each time SGV-caller is started, the software first reads the configuration file "sgv-caller.conf" to obtain the necessary information. Calculations 1-8 in "sgv-caller.conf" correspond to calculation pipelines 1-8, which have different functions as follows:
Pipeline 1 reads the GISAID sequence and metadata files and generates an SGV database from scratch.
Pipeline 2 updates an existing SGV database by comparing the old SGV database with the new GISAID data.
Pipeline 3 reads a FASTA-formatted sequence data and generates an SGV database from scratch.
Pipeline 4 makes the SGV database start from a file containing the raw variations of each sequence.
Pipeline 5 generates the SGV database, which is a subset of the full GISAID data.
Pipeline 6 extract the sequences, which is a subset of the full GISAID FASTA-sequence file.
Pipeline 7 This pipeline extracts the genomic regions of some selected sequences from the GISAID genomes.
Pipeline 8 This pipeline reads a FASTA-formatted protein sequence and maps it to a protein reference sequence.
The GISAID database currently contains approximately 20 million SARS-CoV-2 genomes. Starting the database on a single thread can take more than a month. To speed up the database computation, pipelines 1-3 can be run in parallel with multiple threads, greatly increasing the computation speed.
output_file_name is the run name needed for pipelines 1-8. It should be a string containing only letters, numbers and underscores. No spaces are allowed.
Pipelines 1-5 require a reference genome and its annotation information to annotate nucleotide variations at the codon and amino acid level. By default, SGV-caller uses Wuhan-hu-1 (RefSeq: NC_045512.2) as reference genome. use_default_reference asks if you want to use the default reference genome or not. It takes the value "yes" or "no". If use_default_reference = yes, the default reference genome Wuhan-hu-1 will be used. If use_default_reference = no, the following three files are required to run SGV-caller:
custom_reference_genome
custom_codon_fasta_annotation_file
custom_rna_annotation
custom_reference_genome is the directory and filename of the custom reference genome file. This file should contain a single sequence and its sequence name should not overlap with the sequence name in the fasta data file. custom_codon_fasta_annotation_file is the directory and filename of the custom codon annotation file. Its format is:
S 21563..25384
ORF3a 25393..26220
E 26245..26472
M 26523..27191
nsp10 13025..13441
nsp11 13442..13480
nsp12 13442..13468,13468..16236
…
".." can be used to separate different genetic regions. Note that columns are separated by tabs, NOT spaces. Same for custom_rna_annotation.
custom_rna_annotation is the directory and filename of the custom RNA annotation file. Its format is
5UTR 1..265
3UTR 29675..29903
Non_Coding_ORF1ab_S 21556..21562
Non_Coding_S_ORF3a 25385..25392
…
The calculation value in "sgv-caller.conf" takes values from 1 to 8, corresponding to different pipelines. In the example folder you will find examples and test data of "sgv-caller.conf" for pipelines 1-8.
Pipeline 1: reads the GISAID sequence and metadata files and generates an SGV database from scratch.
calculation = 1
threads = 4
input_GISAID_fasta_genome_file = input_fasta_data_file
input_GISAID_metadata_file = input_metadata_file
Calculation pipeline 1 calculates the local database from the GISAID FASTA sequence data file and metadata file specified by "input_GISAID_fasta_genome_file" and "input_GISAID_metadata_file". Data downloaded from either Genomic Epidemiology or Download Packages can be used directly without unzipping. The GISAID ID, rather than the sequence name, is used as the sequence identifier in this calculation.
Pipeline 2: updates an existing SGV database by comparing the old SGV database with the new GISAID data.
calculation = 2
threads = 4
directory_of_previous_database = directory_of_existing_SGV_database
output_file_name_of_previous_database = output_file_name_of_previous_run
input_new_GISAID_fasta_genome = input_fasta_data_file
input_new_GISAID_metadata = input_metadata_file
Calculation pipeline 2 is used to update an existing database. It takes 4 input values as shown above. It requires the directory of the existing SGV database and its output file name. The GISAID ID, rather than the sequence name, is used as the sequence identifier in this calculation.
calculation = 3
threads = 4
input_fasta_file = multiple_fasta_file
Calculation pipeline 3 is NOT for data downloaded from the GISAID database, but for general sequence data in FASTA format, such as data downloaded from other databases such as GenBank. The sequence names in the multiple_fasta_file are used directly as sequence identifiers. The sequence name should be a string containing only number, character and underscore. No spaces are allowed in sequence names. It is recommended to use the NCBI accession number as the sequence name.
Pipeline 4: makes the SGV database start from a file containing the raw variations of each sequence.
calculation = 4
input_raw_variation_file = input.raw_variants.for_each.all.txt
Calculation pipeline 4 calculates the SGV database from the input.raw_variants.for_each.all.txt file obtained from other runs.
calculation = 5
GISAID_ID_list_file = gisaid_id_list_file
input_GISAID_fasta_genome_file2 = input_fasta_data_file
input_GISAID_metadata_file2 = input_metadata_file
Calculation pipeline 5 reads a "gisaid_id_list_file" and generates a database of GISAID IDs, which is only contained in the "gisaid_id_list_file". It also requires "input_GISAID_fasta_genome_file2" and "input_GISAID_metadata_file2".
calculation = 6
GISAID_ID_list_file = gisaid_id_list_file
input_GISAID_fasta_genome_file2 = input_fasta_data_file
input_GISAID_metadata_file2 = input_metadata_file
Pipeline 6 will extract sequences listed in the gisaid_id_list_file and will NOT generate the SGV database.
calculation = 7
GISAID_ID_list_file = gisaid_id_list_file
input_GISAID_fasta_genome_file2 = input_fasta_data_file
input_GISAID_metadata_file2 = input_metadata_file
gene_annotation = default
#default or the custom gene annotation file name
Pipeline 7 reads a sequence list file gisaid_id_list_file and a gene annotation file gene_annotation. If you specify "default" for gene_annotation, the gene annotation file in the reference folder ("NC_045512.anno.txt") in the software will be used. If you wish to extract the gene regions you have defined, please specify the directory and filename to gene_annotation_file.
calculation = 8
protein_sequence_file = S.pep.fas.gz
protein_reference_sequence_file = NC_045512.2.S.fas
Calculation pipeline 8 reads a protein sequence file as an input file (e.g. S.pep.fas.gz), a protein reference sequence file (e.g. NC_045512.2.S.fas) and extracts the amino acid replacement information directly.
"Output_file_name.raw_variants.for_each.all.txt": the raw variation at the nucleotide level for each genome. "Output_file_name.variations.info.sum.txt": the information for each variation, including SNP locations, reference SNPs, alternative SNPs, reference codons, alternative codons, reference amino acids, alternative amino acids, etc. "Output_file_name.raw_variants.static.txt": the quality of each genome, including the number of differences to the reference genome, the number of undetermined nucleotides of the whole genome, and the number of undetermined nucleotides of the spike protein. "Output_file_name.s.static.txt": the quality of the spike protein region only. "Output_file_name.deletion.txt" and "Output_file_name.insertion.txt": information on deletions and insertions in the genomes analysed. "Output_file_name.count_aa.txt", "Output_file_name.count_codon.txt", "Output_file_name.count_snp.txt": the count of each variation at amino acid, codon and nucleotide level. "Output_file_name.options.txt": the command options that are actually executed during the calculation.
Main output folders:
Output_file_name_genomic_variation
Output_file_name_genomic_variation_ID_unique
Output_file_name_genomic_variation_long_table
Output_file_name_genomic_variation_var_for_each_ID
Each of these folders contains 3 files, aa.txt, codon.txt, snp.txt. These are the summarised genetic variations at amino acid, codon and nucleotide level, with a slightly different data format.
Genetic variation summarised by haplotype. aa.txt in the "Output_file_name_genomic_variation" folder:
nsp2#I120F|nsp12#P323L|S#S477N|S#D614G|N#RG203KR|ORF9c#GE50NE 1 MW155300.1
nsp2#I120F|nsp12#P323L|nsp12#K718N|S#S477N|S#D614G|N#RG203KR|ORF9c#GE50NE 1 MW154115.1
nsp2#I120F|nsp12#P323L|S#S477N|S#D614G|ORF3c#K17R|N#RG203KR|ORF9c#GE50NE 1 MT972245.1
Genetic variation summarised by variation. aa.txt in the "Output_file_name_genomic_variation_ID_unique" folder:
M#L102LZ 1 OU811303.1
N#A211V 1 OD900734.1
N#A220V 7 OU794820.1|OB997169.1|OD944225.1|OD952064.1|OB987694.1|OU030367.1|OU800131.1
N#D288N 1 ON577819.1
N#D343G 1 OW518434.1
contains only two columns, the variation and genomic ID pairs. This format is suitable for analysis by R software or bash commands. Example: aa.txt in "Output_file_name_genomic_variation_long_table" folder:
E#T9I OW513133.1
E#T9I ON556126.1
E#T9I OW464169.1
E#T9I OW506517.1
E#V5F MT994988.1
is the genomic variations for each genome. Example: aa.txt in Output_file_name_genomic_variation_var_for_each_ID folder folder:
FR998558.1 nsp12#P323L|S#D614G|S#Q677H|N#RG203KR|ORF9c#GE50NE
MT451742.1 nsp4#F308Y|ORF3a#G196V|ORF8#L84S|N#P13L|ORF9b#P10S|N#S197L|ORF9c#Q44*
MT810919.1 nsp12#P323L|S#D614G|S#T678I
MT972245.1 nsp2#I120F|nsp12#P323L|S#S477N|S#D614G|ORF3c#K17R|N#RG203KR|ORF9c#GE50NE
In addition to these four folders, SGV-caller also generates three folders summarising the genomic variations for each annotated gene or RNA region, which are summarised in the Output_file_name_aa, Output_file_name_codon and Output_file_name_rna folders.
Example data for each pipeline is prepared in the "examples" folder. An example "sgv-caller.conf" is also prepared with the input data already set for each pipeline. In "sgv-caller.conf" for each example, we have kept only the necessary options for each pipeline. Normally you can leave all the lines in "sgv-caller.conf" as they will not affect the operation of SGV-caller. Any lines in "sgv-caller.conf" that start with "#" are considered annotation lines and will not be read by the software.