deltaRpkm: an R package for a rapid detection of differential gene presence between related bacterial genomes
BMC Bioinformatics volume 20, Article number: 621 (2019)
Comparative genomics has seen the development of many software performing the clustering, polymorphism and gene content analysis of genomes at different phylogenetic levels (isolates, species). These tools rely on de novo assembly and/or multiple alignments that can be computationally intensive for large datasets. With a large number of similar genomes in particular, e.g., in surveillance and outbreak detection, assembling each genome can become a redundant and expensive step in the identification of genes potentially involved in a given clinical feature.
We have developed deltaRpkm, an R package that performs a rapid differential gene presence evaluation between two large groups of closely related genomes. Starting from a standard gene count table, deltaRpkm computes the RPKM per gene per sample, then the inter-group δRPKM values, the corresponding median δRPKM (m) for each gene and the global standard deviation value of m (sm). Genes with m > = 2 ∗ sm (standard deviation s of all the m values) are considered as “differentially present” in the reference genome group. Our simple yet effective method of differential RPKM has been successfully applied in a recent study published by our group (N = 225 genomes of Listeria monocytogenes) (Aguilar-Bultet et al. Front Cell Infect Microbiol 8:20, 2018).
To our knowledge, deltaRpkm is the first tool to propose a straightforward inter-group differential gene presence analysis with large datasets of related genomes, including non-coding genes, and to output directly a list of genes potentially involved in a phenotype.
In comparative genomics the gene presence/absence analysis is commonly performed by multiple alignment calculations on whole genomes or on their subsets as pan-core-genome analysis. Multiple alignment approaches like Mauve  and Mugsy  become quickly very computationally intensive and unsuitable when dealing with increasing number of genomes. For instance, in the case of N = 57 E.coli genomes, Mauve run is not finished after 2 days, while Mugsy needs about 20 h (see ). Pan-core-genome tools like Microscope , Large-Scale Blast Score Ratio (LS-BSR)  require genome assembly and gene prediction steps before performing all-against-all Blast calculations. Roary  performs a clustering of highly similar sequences before executing all-against-all Blast searches only on these subsets of pre-clustered genes, still requiring the assembly and annotation of all genomes . Bacterial Pan-Genome Analysis tool (BPGA)  is fast by clustering the gene sequences like Roary and then aligning them with MUSCLE instead of applying an all-against-all Blast method. Overall, these pan-genome methods run fast on a small scale, e.g., ~ 3 min for BPGA with N = 28 Streptococcus pyogenes samples (genome size ~ 1.8 Mb)  and ~ 6 min for Roary for N = 24 Salmonella enterica, serovar Typhi samples (genome size ~ 4.8 Mb) . However, none of them is practical for larger datasets, e.g., BPGA takes 7 h for 1000 genomes for 4GB of RAM  and Roary produces a pan-genome from 1000 isolates in about 4.5 h, using 13GB of RAM . The above methods are focusing on the protein coding genes, neglecting the non-coding features e.g., small RNA . Other methods like core genome MultiLocus Sequence Typing (cgMLST) are not appropriate for gene presence/absence since the analysis is based on the core-genome, potentially present in all genomes of certain species [9, 10].
Increasing number of studies in human or veterinary clinical genomics, especially those focusing on outbreak detection and tracking, involve a large number of similar genomes to be compared. For such particular cases, we propose a simple yet effective approach using a canonical gene read count table, short-cutting the intensive genome assembly and annotation tasks. Our user-friendly and open-source R package, deltaRpkm, identifies putative genes involved in a given phenotype by inferring their presence/absence from their differential coverage between a reference genome group and a comparison group.
The deltaRpkm pipeline requires as input data metadata and gene read count tables. The read count table can be derived from standard methods like bedtools multicov  based on a reference genome annotation file and the bam files produced by bwa mem . Alternatively, the rapid RNA-seq aligner STAR can be used to obtain the coverage table  (Fig. 1).
Definition of the phenotypic groups
The analysis is centered around a pairwise comparison of gene differential presence between genomes categorized into two different groups according to a selected phenotype: i) a group 1 that shares the phenotype A of the reference genome and ii) a group 2 that does not have the reference phenotype A. This phenotype information per group is provided in the metadata table. The design of the analysis is given in the deltaRpkm::loadMetadata function that loads the grouping criteria of the dataset based on the metadata information.
Conversion of gene read counts to RPKM
The pipeline runs the deltaRpkm::rpkm function to normalize raw read counts with the validated RPKM method (Reads Per Kilobase per Million mapped reads), that takes into account sequencing depth and gene length . For a given sample s of total read counts Ns, the library size correction of read counts (RPMj) corresponds to a scaling factor (scalingFactor) applied to the reads counts per gene (readCountsPerGene), as:
Then, for a given gene j the RPKMj value is computed by weighing in the gene length (geneLength):
Inter-group RPKM values (δRPKM)
For each pairwise comparison of the RPKM values of a gene j between a genome x from group 1 (reference genome) and a genome y from group 2, deltaRpkm::deltarpkm function computes the difference of their RPKM values at gene j (δRPKMj) as:
Selection of genes differentially present in the reference group
The set of genes potentially involved in the selected phenotype correspond to genes that are considered differentially present in the reference genome group, but absent from the comparison group. The deltaRpkm functions to infer those genes are grouped into a main method called deltarpkm::deltaRPKMStats. For each gene j, the median value mj of all its pairwise δRPKM values is calculated, followed by the standard deviation sm of all genes m values. Genes with m > = 2 ∗ sm are considered as present in group 1 of the reference genome and absent from group 2 (Fig. 2). This threshold is relatively stringent and arbitrary, but safer to avoid false positives. Users of deltaRpkm could potentially use the robust Median Absolute Deviation (MAD) as the lower limit to accept a gene differentially present in the reference group. However, this increases the risk of revealing false positives.
Visualisation of the filtered genes
For a more visual evaluation of the selected genes potentially involved in the studied phenotype, deltaRpkm provides a plot function called deltarpkm::rpkmHeatmap which is based on gplots::heatmap.2 method (https://CRAN.R-project.org/package=gplots). This deltaRpkm function plots the RPKM values of the selected genes as a heatmap (Fig. 3). The heatmap color scale is based on the boundaries of the RPKM bimodal distribution (Additional file 1: Figure S1).
The different steps and main functions for a quick start with deltaRpkm are summarized in the Table 1.
The package provides working example datasets of different sizes from Listeria monocytogenes . The complete documentation with more technical details, full tutorial and running R script can be downloaded from the deltaRpkm GitHub project (Fig. 4) and are also provided as Additional files 2 and 3.
The pipeline has been successfully applied in a recent publication  with N = 225 Listeria monocytogenes genomes annotated for their neurovirulence phenotype, as summarized in Fig. 3. Down-sampling tests show the robustness of the method (Additional file 1: Figure S2), with a consistent filtered gene set (Additional file 1: Figure S3). Analyzing a dataset of N = 225 samples takes less than 20 min (Additional file 1: Figure S4) while using less than 4GB of memory (Additional file 1: Figure S5), which makes deltaRpkm an ideal tool for desktop usage. Randomized genome groupings were performed as negative controls, giving shorter and non-robust lists of candidate genes (Additional file 1: Figure S6).
Our strategy in deltaRpkm has two main limitations: 1) the selection and use of a reference strain for read mapping, and consequently the detection of only differential presence of genes in that genome. But this could be overcome by using another strain for the mapping; 2) the non-detection of phenotypic core genes bearing mutations instead of being absent. Direct performance and feature comparisons with other tools are currently difficult, since deltaRpkm is the only one of its kind to perform comparative genomics bypassing the genome assembly and annotation steps. Nevertheless, the Table 2 summarizes the main features of deltaRpkm in comparison to two other nearest tools, BPGA  and Roary .
A powerful feature of deltaRpkm is the inclusion of non-coding genes in contrast to the classical pan-core-genome methods that only target protein-coding genes [4, 6, 7]. The whole genome of the reference is used, and even short non-coding elements are taken into account.
deltaRpkm is a user-friendly R package that makes use of a standard gene counts table to infer a subset of genes potentially involved in a phenotype. The simplicity of its usage, combined with its scalability to large groups of whole genome datasets are the key features of deltaRpkm in the field of comparative genomics.
Availability and requirements
Project name: deltaRpkm.
Project home page: https://github.com/frihaka/deltaRpkm
Operating system(s): Linux, MacOSX, Windows.
Programming language: R.
License: AGPL v3.
Availability of data and materials
The R package deltaRpkm standalone binaries for Linux, MacOS and Windows10 are available are https://github.com/frihaka/deltaRpkm, including tutorial and full documentation.
Reads Per Kilobase per Million mapped reads
Aguilar-Bultet L, Nicholson P, Rychener L, Dreyer M, Gözel B, Origgi FC, et al. Genetic separation of listeria monocytogenes causing central nervous system infections in animals. Front Cell Infect Microbiol. 2018;8:20.
Darling ACE, Mau B, Blattner FR, Perna NT. Mauve: multiple alignment of conserved genomic sequence with rearrangements. Genome Res. 2004;14:1394–403.
Angiuoli SV, Salzberg SL. Mugsy: fast multiple alignment of closely related whole genomes. Bioinformatics. 2011;27:334–42.
Vallenet D, Belda E, Calteau A, Cruveiller S, Engelen S, Lajus A, et al. MicroScope--an integrated microbial resource for the curation and comparative analysis of genomic and metabolic data. Nucleic Acids Res. 2013;41(Database issue):D636–47.
Sahl JW, Caporaso JG, Rasko DA, Keim P. The large-scale blast score ratio (LS-BSR) pipeline: a method to rapidly compare genetic content between bacterial genomes. PeerJ. 2014;2:e332.
Page AJ, Cummins CA, Hunt M, Wong VK, Reuter S, Holden MTG, et al. Roary: rapid large-scale prokaryote pan genome analysis. Bioinformatics. 2015;31:3691–3.
Chaudhari NM, Gupta VK, Dutta C. BPGA- an ultra-fast pan-genome analysis pipeline. Sci Rep. 2016;6:24373.
Vernikos G, Medini D, Riley DR, Tettelin H. Ten years of pan-genome analyses. Curr Opin Microbiol. 2015;23:148–54.
Ruppitsch W, Pietzka A, Prior K, Bletz S, Fernandez HL, Allerberger F, et al. Defining and evaluating a Core genome Multilocus sequence typing scheme for whole-genome sequence-based typing of listeria monocytogenes. J Clin Microbiol. 2015;53:2869–76.
Moura A, Criscuolo A, Pouseele H, Maury MM, Leclercq A, Tarr C, et al. Whole genome-based population biology and epidemiological surveillance of listeria monocytogenes. Nat Microbiol. 2016;2:16185.
Quinlan AR, Hall IM. BEDTools: a flexible suite of utilities for comparing genomic features. Bioinformatics. 2010;26:841–2.
Li H, Durbin R. Fast and accurate long-read alignment with burrows-wheeler transform. Bioinformatics. 2010;26:589–95.
Dobin A, Davis CA, Schlesinger F, Drenkow J, Zaleski C, Jha S, et al. STAR: ultrafast universal RNA-seq aligner. Bioinformatics. 2013;29:15–21.
Mortazavi A, Williams BA, McCue K, Schaeffer L, Wold B. Mapping and quantifying mammalian transcriptomes by RNA-Seq. Nat Methods. 2008;5:621–8.
The authors wish to thank the computing resources of the Interfaculty Bioinformatics Unit of Universities of Fribourg and Bern. The Illumina sequencing was made at the NGS platform of the University of Bern, Switzerland and at GATC Biotech, Konstanz, Germany.
This work was funded by the Swiss National Science Foundation (CRSII3_147692). We declare that the funding agency played no roles in the design of the study and collection, analysis, and interpretation of data and in writing the manuscript.
Ethics approval and consent to participate
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
RPKM values distribution of all genes in the dataset. This can be used to fine tune the heatmap color break parameters. Figure S2. Dataset size effect on the distribution of the δRPKM values. A. Boxplots for datasets from N = 7 to N = 225 samples. The dataset size does not influence the median δRPKM values that are used when computing the differentially present gene selection based on the 2*standard deviation of median δRPKM values. Two datasets are highlighted for illustration, N = 51 samples and N = 225 samples. B. Dataset size effect on threshold value (2*standard deviation) of median δRPKM. Figure S3. The selected differentially present gene set is robust. Downsampling shows that even with small size dataset, the identified genes highly overlap (N = 115) with the datasets of greater size. Figure S4. deltaRpkm performance: dataset size effect on runtime. The whole analysis pipeline with deltaRpkm can be run in less than 20 min in R for a dataset with N = 225 samples of Listeria monocytogenes (~ 3 Mb, ~ 3 K genes). Ubuntu 14.04, R 3.4.4, Intel Core i-4790 CPU @3.60Gzx8. Figure S5. deltaRpkm performance: dataset size effect on memory usage. The whole analysis pipeline with deltaRpkm uses less than 4G of memory in R for a dataset with N = 225 samples of Listeria monocytogenes (~ 3 Mb, ~ 3 K genes). Ubuntu 14.04, R 3.4.4, Intel Core i-4790 CPU @3.60Gzx8. Figure S6. deltaRpkm performance: real (A) versus randomized datasets (B). The gene differential presence gives shorter and non-robust list of genes when using randomized datasets of different sizes. Corrected RPKM.
R usage script for a quick start.
About this article
Cite this article
Akarsu, H., Aguilar-Bultet, L. & Falquet, L. deltaRpkm: an R package for a rapid detection of differential gene presence between related bacterial genomes. BMC Bioinformatics 20, 621 (2019). https://doi.org/10.1186/s12859-019-3234-2
- Comparative genomics
- Differential gene presence/absence