SeroCall can identify and quantitate the capsular serotypes in Illumina whole-genome sequencing samples of S. pneumoniae, calculating abundances of each serotype in mixed cultures. The software is written in Python (compatible with Python 2 or 3), and is freely available under an open source GPLv3 license.
Note (July 30, 2019): The software has just been released, and is still considered in beta testing.
SeroCall has the following dependencies:
- Python 2.7+ or 3+ (tested using 2.7.11 and 3.5.1)
- BWA version >= 0.7.15
- Samtools version >= 1.3 (only used if you change/update the reference database files)
Currently, to install SeroCall, download or clone this github repository, then add the SeroCall directory to your PATH environment variable. Installation through GitHub releases, pip, bioconda and Docker are forthcoming.
usage: serocall [-h] [-o OUTPUT] [-t THREAD] r1file r2file
positional arguments:
r1file R1 FASTQ file (can be gzipped)
r2file R2 FASTQ file (can be gzipped)
optional arguments:
-h, --help show this help message and exit
-o OUTPUT, --output OUTPUT
Prefix for output files [default 'sero']
-t THREAD, --thread THREAD
Number of processors to use [default 1]
In its simplest form, the serocall command just takes the the R1 and R2 fastq files from a paired-end sequencing run. The files can be gzipped or uncompressed, and there are no limitations or patterns to the given filenames.
A "-t #" option can be used to have BWA MEM use multiple threads for processing (BWA MEM is very efficient with parallelizing the alignments, so using multiple threads will reduce the running time of the tool).
SeroCall will generate two output files, a *_counts.txt containing the intermediate "bin" counts, counting how many reads aligned across the regions of the serotype references, and a *_calls.txt containing the final calls. By default, the output file names will use the prefix "sero" (outputting files "sero_counts.txt" and "sero_calls.txt"), but that can be changed using the "-o prefix" option.
The following is an example sero_calls.txt file:
##fileformat=SeroCallv1.0
##NumReads=3420303
##NumUnmapped=129037
##NumGenome=3232094
##NumOther=31298
##NumCapsule=27874
#SEROTYPE PERCENTAGE
19F 61.8%
01 38.2%
The first line of the output is a file format tag. The other lines, before the "#SEROTYPE..." header, contain metrics about the read data, including the following:
- NumReads - the number of input reads
- NumUnmapped - number of reads not mapped to serotype or S. pneumoniae genomes
- NumGenome - number of reads mapped to the S. pneumoniae genomes (not the serotype sequences)
- NumOther - number of other alignments, including partial, high error and chimeric alignments
- NumCapsule - number of reads well-aligned to the serotype capsular sequences (i.e., used to call serotypes)
The "#SEROTYPE..." header and following lines are tab-delimited lines reporting the serotypes called and their percentages. If no lines appear after the header, no serotypes were called from the data.
The sero_counts.txt file is technically an intermediate file, but can be useful for post-analysis QC of the calls. The initial lines of an example file are the following:
##fileformat=BinCountsv1.0
##NumReads=3420303
##NumUnmapped=129037
##NumGenome=3232094
##NumOther=31298
##NumCapsule=27874
#SAMPLE SEROTYPE START END TOTAL UNIQUE
sero 01 1 500 373 369
sero 01 501 1000 546 508
sero 01 1001 1500 612 584
sero 01 1501 2000 640 339
sero 01 2001 2500 512 322
sero 01 2501 3000 556 556
sero 01 3001 3500 568 551
sero 01 3501 4000 669 669
sero 01 4001 4500 637 637
sero 01 4501 5000 852 852
sero 01 5001 5500 775 775
sero 01 5501 6000 698 698
sero 01 6001 6500 655 655
sero 01 6501 7000 729 729
sero 01 7001 7500 786 786
sero 01 7501 8000 732 732
sero 01 8001 8500 763 763
sero 01 8501 9000 780 780
sero 01 9001 9500 729 729
.
.
.
After the header lines, the rest of the file is a tab-delimited file identifying the regions of the serotype sequences, and then reporting the total counts of aligned reads, and the count of uniquely aligning reads in each region. The final lines in the file contain an extra column with a "diff" identifier, marking them as difference locations for distinguishing closely related serogroups.
WARNING: This functionality has not been tested externally, does require some (to a lot of) care when making changes, and should be needed by only a small subset of users.
The current software has pre-built databases containing the same serotype sequences as the PneumoCaT v1.2 database, and building the database is not required in order to use the serocall command. The package does provide the ability to change and update the references used by SeroCall, but this has not been extensively tested.
The builddata sub-directory contains the data files and scripts necessary to rebuild the SeroCall database files (found in the "data" sub-directory of the repository). In that sub-directory, the serotype capsular sequences are found in the "serotypes" sub-directory, the S. pneumoniae references are found in the "genomes" sub-directory and the difference locations are found in the "pcat_diffs.txt" file.
To add a new serotype capsular sequence to the database, create a FASTA file containing only that sequence, and then place it in the serotypes sub-directory. IMPORTANT: The file must end with ".fa" and the name of the file must match the accession given on the FASTA header. So, for a new serotype like "10X", the name of the file should be "10X.fa" and the first line of the file should be ">10X ..." where after the first space on the line, any text is allowed.
To add a new genome reference to the database, create a masked FASTA file of the reference sequence and place it in the genomes sub-directory. Any name can be given to the file, however the accession on the FASTA header must begin with "Spneumo_". IMPORTANT: The capsular sequence must be masked from the genome sequence (or all of the reads will align to it, confusing SeroCall). To do that, megablast/blat/minimap2 the serotype sequences against the genome reference, identify any region of the genome with very high identity (>98%) longer than 100 bp, and replace those regions' bases with N's in the FASTA file. (Remember, the goal is not to match to this sequence, but to use it as a decoy so that genomic reads don't align to the serotype sequences.)
If you believe you want to make changes/additions to the difference locations, please email [email protected] for instructions on how to do that (there are a number of subtle issues that must be accounted for, in order for the software to properly identify serotypes using these difference).
Once the changes have been made to these files, then the command "bash builddata.sh" command will rebuild the database files and install them into the top-level "data" directory, for use by the serocall command. Things to note when running this command:
- Both bwa and samtools are required for this step, so they should be in your PATH.
- All of the initial construction of the files occurs in the "build" directory, then a final transfer occurs into the top-level "data" directory.
- A backup of the existing top-level "data" directory files is made, placed in "data/bak". (But, only a single backup is made, so if builddata.sh is run multiple times, the previous backup is removed.)
- This command will download a pair of FASTQ files, totalling 1.5 GB, and store them in the "genomeOnly" sub-directory the first time it is run. These files contain a real run of a capsule-free sample, and is used to identify those serotype regions where genomic reads may align into the serotypes.
To cite the SeroCall package, please use the following citation:
Knight J, Dunne E, Mulholland E, Saha S, Satzke C, Tothpal A, Weinberger D. Determining the Serotype Composition of Mixed Samples of Pneumococcus using Whole Genome Sequencing. To be submitted.