BingleSeq: a user-friendly R package for bulk and single-cell RNA-Seq data analysis

Implementation

BingleSeq is based on shiny (Chang et al., 2016) and it is composed of a multi-tabbed UI, built as separate shiny modules with efficiency, code readability and reusability in mind. Each module (tab) corresponds to a key step in the typical Bulk and scRNA-Seq analysis pipelines (Fig. 1). Modules are generated only upon reaching a given step of the analyses which ensures efficiency and speed despite the complexity of the application. BingleSeq’s UI components (e.g., plots, tables, tabs, etc.) make use of shiny’s “reactivity” property and these components are automatically updated upon user input or any change related to a given component. The development of BingleSeq focused on providing a flexible and intuitive user experience. As such, BingleSeq implements three state-of-the art package options for each analysis, while it also enables users to tailor key analysis parameters according to their experiment.

User experience is complemented with customizable and interactive tables and plots as well as analysis-related tips and pop-up messages to guide the correct execution of the pipelines. Plot interactivity was achieved via conversion from ggplot to ggplotly (Sievert, 2020).

Bulk RNA-seq steps and features

To begin the DE analysis of Bulk RNA-Seq data, a count table and a metadata table must be loaded in the appropriate formats (Fig. S1). Genes can then be filtered according to counts per million (CPM), Max, or Median thresholds. Furthermore, batch effect correction can be performed with Harman and ComBat packages (Leek et al., 2012; Oytam et al., 2016) (Fig. S2).

Subsequent to quality control, users can investigate differentially expressed genes (DEGs) using three state-of-the-art packages: DESeq2 (Love, Huber & Anders, 2014), edgeR (Robinson, McCarthy & Smyth, 2010), and limma (Ritchie et al., 2015).

Upon obtaining DE results, users can visualize them using interactive plotting techniques (Figs. 2A–2E). These include: a PCA plot used to provide insights about the relationship between samples; Barchart plot supplemented with a summary table used to summarize the up- and downregulated DEGs; Volcano and MA plots that assess the relationship between fold change (FC) versus significance and average expression; and Heatmaps which are arguably the most versatile and informative type of visualization technique when looking at DE results.

BingleSeq’s visualization techniques were implemented with customization in mind and users can specify parameters such as p-value threshold and fold-change threshold, among others. Due to their versatility, heatmaps were designed as BingleSeq’s most customizable plotting component (Fig. 2F).

To assess BingleSeq’s Bulk RNA-Seq pipeline, we used a synthetic dataset generated with the compcodeR package (Soneson, 2014). We also tested it on a real data set looking at DEGs between HSV-1 infected control and interferon B treatment (McFarlane et al., 2019)—see Article S1.

scRNA-seq pipeline steps and features

BingleSeq’s scRNA-Seq part is based on Seurat’s scRNA-Seq pipeline and visualizations (Satija et al., 2015). Nonetheless, clustering can also be performed with monocle (Trapnell et al., 2014) and SC3 (Kiselev et al., 2017) packages.

To begin scRNA-Seq analysis, data can be supplied in two formats—“Cell Ranger 10× Genomics” data and a count table in a predefined format (Fig. S3). Once the data is loaded, users can filter unwanted cells and features (Fig. S4). The next step is to normalize the data and BingleSeq provides two Seurat-supplied normalization methods—“LogNormalize” and “Relative counts”. Simultaneously with normalization and scaling, the highly variable features within the dataset are identified and selected for clustering with Seurat as a way to minimize noise. The highly variable features are visualized using an interactive scatter plot (Fig. S4).

Following normalization, the “Clustering” tab is generated (Fig. 3) which provides a high degree of control over the different steps of the analysis (Fig. 3A). Additionally, general tips are provided, such as clustering advice provided for each package (Fig. 3B). Moreover, each of the visualization techniques in this tab are interactive to further facilitate the interpretation of results. First, the clustering tab requires the generation of pre-clustering prerequisites such as scaling the data and dimensionality reduction with Principal Component Analysis (PCA). Then an elbow plot is returned which is used to determine the dimensionality of the dataset (Fig. 3C)—this is essential for excluding noise when clustering with Seurat and monocle. PC heatmaps (Fig. 3D) are also available as a further tool for PC Selection. Once the data is filtered and transformed, users can proceed to unsupervised clustering with Seurat, SC3, and monocle. The primary way to visualize clustering results is via t-distributed stochastic neighbor embedding (tSNE) plots (Figs. 3E–3G)—a method designed for the purpose of visualizing high dimensional datasets (Van der Maaten & Hinton, 2008).

Following clustering, DE analysis can be conducted using Seurat’s inbuilt testing methods to identify marker genes. The implemented Seurat DE methods include: Student’s T test, Wilcoxon Rank Sum test, DESEq2 (Love, Huber & Anders, 2014), and MAST package (Finak et al., 2015). DE results and specific marker genes can then be visualized using Seurat’s inbuilt plots. These plots include a cluster heatmap and visualizations for the exploration of specific genes via Violin, Feature, and Ridge plots (Fig. 4).

An evaluation of BingleSeq’s scRNA-Seq pipeline was performed by reproducing and extending the results of Seurat’s online tutorial (https://satijalab.org/seurat/v3.0/pbmc3k_tutorial.html)—see Article S2. The tutorial is based on a 10× Genomics dataset of 2,700 Peripheral Blood Mononuclear Cells (PBMCs) with ~69,000 reads per cell. The data set is available at https://support.10xgenomics.com/single-cell-gene-expression/datasets/1.1.0/pbmc3k.

Functional annotation

For both Bulk and scRNA-Seq pipelines, BingleSeq enables the functional annotation of results using GOseq (Young et al., 2010) and transcription factor (TF) and pathway activity inference tools DoRothEA (Garcia-Alonso et al., 2019) and PROGENy (Schubert et al., 2018; Holland, Szalai & Saez-Rodriguez, 2019; Holland et al., 2020).

GOSeq is implemented in the “Gene Ontology” tab and it enables users to obtain results from KEGG pathway analysis and three types of GO categories, including “Cellular Component”, “Molecular Function”, and “Biological Function” (Fig. S5). The “Gene Ontology” tab can also be used to generate top 10 GO term histograms and to obtain additional information about a given GO term using the “GO.db” package (Carlson et al., 2019) (Fig. S5). Note that BingleSeq supports human, mouse, drosophila, zebrafish, and Escherichia coli K12 strain genomes (Carlson, 2019a, 2019b, 2019c, 2019d, 2019e) for analysis with GOSeq.

In BingleSeq, the signed TF-target interactions from DoRothEA are coupled to the analytic Rank-based Enrichment Analysis (aREA) method from the viper package (Alvarez et al., 2016) and together they enable TF activity estimation (Figs. S6A and S6B). On the other hand, PROGENy uses a linear model to estimate the activity of 14 signaling pathways, built on downstream expression patterns from a compendium of perturbation experiments (Figs. S6C and S6D). DoRothEA and PROGENy are footprint-based tools and they infer activity from the expression of molecules considered to be downstream of a given TF or pathway, respectively. Both of these tools are implemented in the “Footprint Analysis” tab and are available for human and mouse data.

DE package comparison and rank-based consensus

BingleSeq supplies an option to assess the agreement between the different DE analysis methods in the form of a Venn diagram (Fig. S7A). Moreover, BingleSeq provides a rank-based consensus approach to improve confidence in DE results (Costa-Silva, Domingues & Lopes, 2017; Guo et al., 2014; Moulos & Hatzis, 2015)—see Fig. S7B. In the context of the Bulk RNA-Seq pipeline the overlap in DEGs is assessed on the results obtained using DESeq2, edgeR, and limma packages. In the case of scRNA-Seq, this is done using three of Seurat’s inbuilt DE methods—MAST, Wilcoxon Rank Sum Test, and Student’s T test.

Inbuilt bulk and single cell RNA-seq example datasets

BingleSeq also provides inbuilt Bulk RNA-Seq and scRNA-Seq test data. Bulk RNA-Seq data is represented by a 3-sample contrast between HSV-1 infected control and interferon B treatment (McFarlane et al., 2019). The data files are available from European Nucleotide Archive (ENA) under accession number PRJEB27501. The single-cell RNA-Seq example is a Cell Ranger 10× Genomics dataset looking at filtered data of 2,700 human peripheral blood mononuclear cells.