Srst2

From wiki
Revision as of 16:13, 19 July 2016 by Rf (talk | contribs)
Jump to: navigation, search

Introduction

This is a python based tool with two main dependencies (samtools and bowtie2) which 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.

Command lines

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

Before anything else, you will need to get hold of an MLST scheme. This will be the fasta sequences of all the alleles of each housekeeping gene which has been slecte for the species in question. In the example, it's obtained as follows:

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

It's unclear what the "#1" means but it probably refers to certain strain of Ecoli.

If looking for Staphylococcus Aureus, the command would be

getmlst.py --species "Staphylococcus Aureus"

For Ecoli, you will get 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. In srst2 terms, this file is the "MLST database".
  • ecoli.txt, a listing of predetermined types, i.e. the sequence types, together with the corresponding alleles of each of the housekeeping genes they should match up with, if they are to be considered "belonging to that type". In srst2 terms these are the "MLST definitions" which we'll come across later.
  • a log file more or less verifying these downloads and their URLs.

In essence therefore, we get the different allele sequences of the housekeeping genes as a "fasta" suffixed file, and then a "txt"-suffixed file with the specifications of each ST in terms of the housekeeping alleles. The allele number is concatentated to the housekeeping gene id line, and often a delimiter must be defined as to how it is concatenated.

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, a file of fast sequences for each of the variants of the housekeeping gene. The number of varaitns for each is unsurprisingly different. adk7 has some 580, while fumC has over 800 variants.
  • --mlst_definitions, defines a file give all the types in terms of their combination of variant number 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 this example, we specify three pairs of fastq files using wild cards and with the --input_pe option.
  • A separate srsts2 instance will be launched within an internal geenrated qsub jobscript for each pair of fastq files.
  • This will create, in this example, three separate job numbers whihc should be visible using the 'qstat command.
  • 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.