deltaRpkm: an R package for a rapid detection of differential gene presence between related bacterial genomes

Background 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. Results 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). Conclusions 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.


Background
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 [2] and Mugsy [3] 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 [3]). Pan-core-genome tools like Microscope [4], Large-Scale Blast Score Ratio (LS-BSR) [5] require genome assembly and gene prediction steps before performing all-against-all Blast calculations. Roary [6] performs a clustering of highly similar sequences before executing allagainst-all Blast searches only on these subsets of preclustered genes, still requiring the assembly and annotation of all genomes [6]. Bacterial Pan-Genome Analysis tool (BPGA) [7] 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 pangenome methods run fast on a small scale, e.g.,~3 min for BPGA with N = 28 Streptococcus pyogenes samples (genome size~1.8 Mb) [7] and~6 min for Roary for N = 24 Salmonella enterica, serovar Typhi samples (genome size~4.8 Mb) [6]. However, none of them is practical for larger datasets, e.g., BPGA takes 7 h for 1000 genomes for 4GB of RAM [7] and Roary produces a pan-genome from 1000 isolates in about 4.5 h, using 13GB of RAM [6]. The above methods are focusing on the protein coding genes, neglecting the non-coding features e.g., small RNA [8]. 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.

Input files
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 [11] based on a reference genome annotation file and the bam files produced by bwa mem [12]. Alternatively, the rapid RNA-seq aligner STAR can be used to obtain the coverage table [13] (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 del-taRpkm::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 [14]. 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: Fig. 1 Overview of a deltaRpkm workflow. Black arrows indicate the main pipeline; dotted arrows show an alternative route with STAR. The package is written in R and takes as input a canonical coverage table, plus the design information given by the user as a metadata table. The strength of deltaRpkm relies on bypassing the tedious assembly and annotation steps typical of comparative genomics. Instead, deltaRpkm uses a basic gene read counts table (based on the mapping against a reference genome) to compute inter-group differential RPKM values per gene and outputs a list of candidate genes as present in the samples of the reference genome group (and absent from the comparison group) scalingFactor ¼ N s 10 6

RPM j ¼ readsCountsPerGene scalingFactor
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 m j of all its pairwise δRPKM values is calculated, followed by the standard deviation s m of all genes m values. Genes with m > = 2 * s m 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.

Tutorial
The package provides working example datasets of different sizes from Listeria monocytogenes [1]. 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.

Results
The pipeline has been successfully applied in a recent publication [1] 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).

Discussion
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 Fig. 3 Heatmap of the RPKM distribution of the selected genes. These genes are considered as differentially present between group 1 (samples that have the same phenotype as the reference genome) and group 2 of samples. A dataset of N = 51 of Listeria monocytogenes genomes is represented in this figure Table 1 Main functions for a differential gene presence/absence analysis with deltaRpkm. Functions are listed in the chronological order of usage  Table 2 summarizes the main features of deltaRpkm in comparison to two other nearest tools, BPGA [7] and Roary [6]. A powerful feature of deltaRpkm is the inclusion of non-coding genes in contrast to the classical pan-coregenome 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.
Conclusions 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 del-taRpkm in the field of comparative genomics.
Additional file 1: Figure S1. RPKM values distribution of all genes in the dataset. This can be used to fine tune the heatmap color break parameters. Figure S2.  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.
Additional file 2. R usage script for a quick start.