Srst2

Introduction

This python based tool with two main dependencies (samtools and bowtie2) carries out mapping of short reads to detect three overall targets:

1. genes
2. alleles
3. multi-locus sequence types (MLST)

from WGS data (which we can take to be NGS short reads).

A smaller number of loci, i.e. 7, are used to divide up the population.

These loci come in the shape of entire housekeeping genes that all the species and isolates are bound to have.

This proves to be robust unlike the alternative of assembling the genomes de-novo, which, when dealing with 100 or 1000 bacterial genomes (a typical workload in bacteria), is a major issue.

Example command lines

(the following uses the example procedure which is described on the project's github page).

In the beginning, you will need to get hold of an MLST scheme. In the example, it's obtained as follows:

getmlst.py --species "Escherichia coli#1"

This will get you the following files:

• The allele sequences of the 7 housekeeping genes: sadk.tfa, fumC.tfa, gyrB.tfa, icd.tfa, mdh.tfa, purA.tfa, recA.tfa.
• Escherichia_coli#1.fasta, seemingly a concatenated version of the above. In fact the above are not used, so this file must represent the above.
• ecoli.txt, a listing of STs together with their allelic profile for each of the 7 housekeeping genes.
• a log file more or less verifying these downloads and their URLs.

Each housekeeping gene has about 200-300 sequences (i.e. allele-sequences) of length in the order of roughly about 400bp.

Then we need the short read pairs. Each pair probably represents a sample. The example uses five pairs. It's unclear what each represents. First off, a low coverage pair is chosen which will not give srst2 much signal:

srst2 --input_pe ERR024070*.fastq.gz --output shigella1 --log --save_scores --mlst_db Escherichia_coli#1.fasta --mlst_definitions ecoli.txt --gene_db ARGannot.fasta

Explanation:

• Up to --save_scores the command line clearly takes in input fastq files, and defines an output directory. The final three options are now explained.
• --mlst_db, this is for specifying the MLST database, usually a fasta file with the allele sequences for certain marker sequences. In srst2 this will be the concatenated fasta of different versions of the 7 housekeeping genes which correspond to the different allele sequences of the gene.

• --mlst_definitions, defines a file give all the types in terms of their combination of the seven housekeeping genes.

• --gene_db, this is a fasta file of genes that confer antibiotic resistance. It allows for gene-typing.

Another important option, not invoked here, is

• --mlst_delimiter, the default (as suited the example) is '-', but you can set it to some other character, remembering the single quotes.

A somewhat diverting option is

• -output which actually is for giving the the prefix name for the various files that will be output.

The end of the output from this command will give a idea as to how rich the results, and they will be poor if there are few read pairs that aligned. Thsi will be true in the case of the example's ERR024070 read pair, which was very low coverage. A good deal of files are also output from this stage:

• .bt2 and .fai files: these are from the bowtie2 indexing and aligning stage which occurs on the reads against both databases.
• .sorted.bam and .pileup files for both alignment operations, which are the results of the alignment operated on by samtools.
• Two .scores files corresponding to the above.
• Three result files, named appropriately with .txt extensions. Only one comes from the MLST database while the gene database gives two, one is a summary simply and the other gives full details. The subname "fullgenes" does not refer to full or partial genes but rather the report itself, whether it is a summary or not. The summary file will add columns for each of ARGs identified in the read-pair.
• a .log file detailing these operations.

The example's main concern is the STs that have been called. The file for this result has name of template [outputprefix]__mlst__[db]__results.txt. While a question mark means low certainty, an NF means not found. This is what we have in this case, an NF. ND would have meant no definition.

Having tried out the first read pair, the procedure continues with the second read pair building upon previous results, in this case the low coverage pair just analysed. Here is the command:

srst2 --input_pe ERR028678*.fastq.gz --output shigella2 --log --mlst_db Escherichia_coli#1.fasta --mlst_definitions ecoli.txt --gene_db ARGannot.fasta --prev_output shigella1__compiledResults.txt

Notes:

• --prev_output, here you can see the option that allows usage of a previous result. This is sort of pooling of the read-pairs, it's unclear what the experimental conditions could be for this.

So two read-pairs have been processed. Multi-pair processing can also be done with a single command-line. Though clearly they will be done independently, without prior consideration unlike the first two pairs, where the second pair built upon the results of the first. However, the results can be compiled afterwards as will be shown.

The multi-pair processing command-line is as follows:

srst2 --input_pe ERR024082*.fastq.gz ERR028690*.fastq.gz ERR024619*.fastq.gz --output shigella3 --log --mlst_db Escherichia_coli#1.fasta --mlst_definitions ecoli.txt --gene_db ARGannot.fasta

Note how the output prefix is shigella3, our third processing command-line.

Then the compilation step is as follows:

srst2 --output all --prev_output shigella2__compiledResults.txt shigella3__compiledResults.txt

Command-line notes:

• -prev_output here accepts several previous result files
• -output again gives the prefix name, so that output file will be called all__compiledResults.txt.

Parallel srst2

The srst2 team put together a script called slurm_srst2 for running on a cluster administered by the queue manager SLURM. Marvin uses a different queue manager called Gridengine (qsub) but this script has been modified so that the Gridengine queue manager will accept it and therefore run the datasets in parallel (essentially, each job becomes a separate job, with jobs launching simultaneously to achieve parallelism. This script is called qsub_srst2 and is similar, but with a few differences, to the srst2 command. A typical command line is as follows:

 qsub_srst2.py --input_pe ERR024082*.fastq.gz ERR028690*.fastq.gz ERR024619*.fastq.gz --output shi_para_0 --queue lowmemory.q --other_args "--mlst_db Escherichia_coli#1.fasta --mlst_definitions ecoli.txt --gene_db ARGannot.fasta"

Explanation

• in the --input_pe we specify three pairs of input files which will be operated upon in three separate batches, in parallel. It will create three separate job numbers.
• the preceding steps of bowtie and samtools steps will not be done in parallel
• the --other_args (note the underscore) is a way of feeding the extra arguments that the normal srst2 script used. Note also how these must be presented in quotations marks, as a single string.

Glossary

• ST: sequence type ... can be seen as a category a certain sequence can be said to belong to.
• MLST: MultiLocus Sequence Type
• ARG: Antibiotic Resistance Gene, a gene that imbues the microbe with resistance.