Development of a bioinformatics platform for analysis of quantitative transcriptomics and proteomics data: the OMnalysis

Software packages and implementation

OMnalysis is an interactive R Shiny based web application composed of multiple sectioned user interface (UI) in the form of tabs. It is built to perform exploration of differential expression data efficiently and iteratively. R packages used for the development of the UI and its components are as follows: R Shiny version 1.6.0 (Chang et al., 2017), flexdashboard version 0.5.2 (Iannone, Allaire & Borges, 2018), Shiny Themes version 1.2.0 (Chang et al., 2018), rmarkdown version 2.8 (Allaire et al., 2020), knitr version 1.33 (Xie, 2019) and Shiny dashboard version 0.7.1 (Chang, Ribeiro & Barbara, 2019). Each sectioned UI is further divided into interactive and display panels. The interactive panel works on the Shiny’s reactivity property, which automatically updates the values in the output panel when the user interacts or changes the input components (plots, tables, actions, etc.). The biological analysis is supported by the following R packages: biomaRt version 2.46.3 (Durinck et al., 2009), clusterProfiler version 3.18.1 (Yu et al., 2012), reactomePA version 1.34.0 (Yu & He, 2016), reactome.db version 1.74.0 (Ligtenberg, 2019), pathview version 1.30.1(Luo & Brouwer, 2013), SPIA version 2.42.0 (Tarca et al., 2009), SBGNview version 1.4.1 (Dong et al., 2021), STRINGdb version 2.2.2 (Szklarczyk et al., 2019), org.Hs.eg.db, org.Gg.eg.db, org.Ss.eg.db, org.Bt.eg.db, (org.Mm.eg.db), (org.Rn.eg.db), (org.Cf.eg.db), (org.Dm.eg.db) and (org.Ce.eg.db)version 3.12.0 (Carlson, 2019). Visualization of the results from the analysis is backed by the following packages: EnhancedVolcano version 1.8.0 (Blighe, Rana & Lewis, 2020), gplots version 3.1.1 (Warnes et al., 2016), ggbiplot version 0.55 (Vu, 2016), ggplot2 version 3.3.3 (Wickham, Chang & Wickham, 2016), VennDiagram version 1.6.20 (Chen & Boutros, 2011), wordcloud version 2.6 (Fellows, 2018), dplyr version 1.0.5 (Wickham et al., 2015) and DT version 0.18 (Xie, Cheng & Tan, 2018).

To enhance the user experience and proper execution of the OMnalysis pipeline (Fig. 1), peer-reviewed R packages were streamlined in multi-tabbed UI as follows: Upload data, PCA, Plots, Statistical filtering, GO enrichment analysis, GO heatmaps, Pathway enrichment analysis, Enriched pathway visualization, Literature info and Help.

The web application is designed to analyze two types of quantitative omics, the transcriptomics and proteomics. For transcriptomics, RNA-Seq data generated previously by us and deposited in ArrayExpress was used (http://www.ebi.ac.uk/arrayexpress). Expression analysis was performed on human brain microvascular endothelial cells (hBMEC) induced with various pathogens: Borrelia burgdorferi (Treatment1, Table S1, retrieved from ArrayExpress accession number E-MTAB-8053), Neisseria meningitidis (Treatment2, Table S2, E-MTAB-8008), Streptococcus pneumoniae (Treatment3, Table S3, E-MTAB-8054), and West Nile Virus (Treatment4, Table S4, E-MTAB-8052). Table S1 to Table S4 are in text format processed from the TSV file generated from the edgeR’s glmTreat function (Robinson, McCarthy & Smyth, 2010). Three columns (logFC, logCPM, and P value) from each of those supplemental files were copied to make a master file in CSV format (Table S5). For the sake of simplicity, Table S5 is explained in Table 1.

In the case of proteomics, the data matrix was generated from one of the differential abundance analysis software Perseus version 1.6.15 (Tyanova et al., 2016). To check the functionality, we retrieved the .xlsx format file from the experiment performed to quantify protein abundance in the milk whey collected at different time points from the cow with Streptococcus uberis infection (Mudaliar et al., 2016). The columns in this data matrix were arranged in the following order: UniProt ID, FDR-adjusted P-value, and Fold Change in an excel file for each experimental condition (4 time points in this case, designated as Treatment1, Treatment2, Treatment3, and Treatment4; Table S6).

Data modification and ID conversion

We used the read.csv function of Utils package version 3.6.2 (Team, 2013) to upload the CSV format file (Table S5) to the OMnalysis using a tab “differentially expressed example data”. Whereas, for proteomics, we used the import_list function of rio package version 0.5.26 (Chan et al., 2018) to upload the data (Table S6) using a tab “proteomics abundance example data”. For proteomics data, three functions were used to convert data from Table S6 to make input table for the OMnalysis. First, the rio package was used to convert Treatment sheets to Treatment column. Second, the boldlog2 function of base R was used to transform the Fold Change column values to logFC (log Fold Change), and the third, Colnames function was used to change the treatment column name to Treatments. The duplicate proteins in the treatments were identified using the group_by function of dplyr version 1.0.5 (Wickham et al., 2015) and the mean function of base R to obtain the mean of their log fold change value. Such conversion is not necessary in case of transcriptomic data.

Transcriptomic data matrix contains Ensembl IDs, while proteomic data comes with UniProt IDs. To convert these IDs into five different ID types (Ensemble gene ID, gene name, HGNC symbol, gene description, and UniProtKB/Swiss-Prot ID) we used the getBM function of biomaRt package version 2.46.3 (Durinck et al., 2009) to fetch the latest information from the Ensembl database (Yates et al., 2020). We have incorporated the possibility of ID conversion for 9 species (Human, Chicken, Pig, Cow, Mouse, Rat, Dog, Drosophilla melanogaster and C.elegans) in OMnalysis.

Principle component analysis (PCA)

To perform PCA, following the data matrix upload, the UI PCA tab (Fig. 2A) was included. In the interactive panel (Fig. 2B) for PCA, checkboxes and dropdown menu were inserted for Variable PCA and Biplot PCA. The fviz_pca_var function of factoextra package version 1.0.7 (Kassambara & Mundt, 2017) and prcomp function of stat version 3.6.2 (Team, 2013) were used for Variable PCA and Biplot PCA, respectively. The biplot function of ggbiplot version 0.55 (Vu, 2016) was used for the visualization (example data is shown in Figs. 2C and 2D).

Plots

We have created UI tab Plots (Fig. 3A) for generation of the scatter and volcano plots. An interactive panel (Fig. 3B) was created using flexdashboard version 0.5.2 (Iannone, Allaire & Borges, 2018), which accommodates the following options to users: selection checkboxes to compare the treatments and numeric input boxes to add values for P value and log fold changes (P value and FC-cutoff in Fig. 3B). We also introduced three checkboxes in the interactive panel for Scatter and Volcano plot for transcriptomics, whereas, Volcano plot checkbox for proteomics. Multiple image format drop-down menu and dimension correction numeric input options were also included to customize the generated plot (Fig. 3B). EnhancedVolcano package version 1.8.0 (Blighe, Rana & Lewis, 2020) and ggplot2 package version 3.3.3 (Wickham, Chang & Wickham, 2016) were used to generate an example volcano plot (Fig. 3C) and scatter plot (Fig. 3D) from DEGs or DEPs, respectively.

Statistical filtering

Following the PCA and plots, we included the Statistical filtering option in the UI tab (Fig. 4A). The interactive panel in this UI was populated with checkboxes to select treatments (where user can select the treatments for filtering and comparison). A dropdown menu “Omics Type” was added to select transcriptomics or proteomics. A numeric input box “Statistical filtering” was added to insert values for cutoff. Various checkboxes under the Venn Diagram and Histogram were added to plot the graphs based on cutoff values (Fig. 4B). The VennDiagram package version 1.6.20 (Chen & Boutros, 2011) and ggplot2 package version 3.3.3 (Wickham, Chang & Wickham, 2016) were used to plot Venn diagram (Figs. 4C and 4D) and histogram (Figs. 4E and 4F), respectively.

Gene ontology analysis

Once the data matrix was statistically filtered, we used clusterProfiler version 3.18.1 of Bioconductor packages (Yu et al., 2012) to obtain the functional interpretation of the significantly expressed genes or proteins. Two enrichment analysis functions of the clusterProfiler were used, the first, EnrichGO function on genes or proteins to perform over representation analysis (ORA) and the second gseGO function on sorted genes or proteins with respect to logFC values to perform gene set enrichment analysis (GSEA) (Fig. S1). To support the enrichment analysis, AnnotationDbi version 1.52.0 of Bioconductor databases (Pagès et al., 2020) was used. Mark Carlson species-specific genome-wide annotation databases version 3.12.0 was used for human (org.Hs.eg.db), chicken (org.Gg.eg.db), pig (org.Ss.eg.db), and cattle (org.Bt.eg.db), mouse (org.Mm.eg.db), rat (org.Rn.eg.db), dog (org.Cf.eg.db), Drosophila melanogaster (org.Dm.eg.db) and C.elegans (org.Ce.eg.db) (Carlson, 2019). P-value cutoff input and PAdjust function of ClusterProfiler with seven multiple hypotheses testing correction methods ( Fig. S1) were used to avoid the influence of false-positive results on the overall enrichment analysis.

Heatmap and word cloud

In the GO heatmaps UI tab (Fig. 5A) we included an interactive panel (Fig. 5B), which was populated with checkboxes for selection of heatmap visualization method, treatments checkboxes, numerical inputs boxes to adjust the font and color key size, word cloud checkbox to plot a word cloud, etc. The wordcloud package version 2.6 (Fellows, 2018) was employed to visualize a word cloud (Fig. 5C), while heatmap.2 function of gplots version 3.1.1 (Warnes et al., 2016) and rainbow function of the R base version 4.0.3 (RStudioTeam, 2015) were used to generate an example heatmap (Fig. 5D).

Pathway analysis

To get the mechanistic insight from the list of DEGs or DEPs, four different pathway analysis methods were explored (Fig. S2). To enrich the DEGs or DEPs against KEGG pathway database (Kanehisa & Goto, 2000), two functions namely enrichKEGG and gseKEGG from cluster profiler’s version 3.18.1 (Yu et al., 2012) were used.

The network topology analysis (NTA) (Alexeyenko et al., 2012), was used against reference databases such as biocarta (Rouillard et al., 2016), panther (Thomas et al., 2003), NCI- nature pathway interaction database (Anthony et al., 2011), and pharmGKB (Klein & Altman, 2004). For NTA, we used three function in R base version 4.0.3 namely merge, cbind and gsub (RStudioTeam, 2015) to arrange the data matrix of DEGs or DEPs. Then, we used the graphite pathway function of graphite package version 1.36.0 (GRAPH Interaction from pathway Topological Environment) (Sales et al., 2012) to prepare the reference pathway database. The graphite’s runSPIA function was used to perform network topology analysis using the four reference databases mentioned above.

An enrichpathway function of ReactomePA version 1.34.0 (Yu & He, 2016) was used to perform pathway analysis using the Reactome pathway database version 1.74.0 (Croft et al., 2011).

A $new function in STRINGdb version 2.2.2 (Szklarczyk et al., 2019) was used to assign the species, score threshold, and input directory. We used the stringdb$map function of the STRINGdb package to map the DEGs or DEPs against several databases (GO annotation, KEGG pathways, PubMed publications, Pfam domains, InterPro domains, UniProt Keywords SMART domains).

Enriched pathway visualization

Enriched pathway visualization tab of UI (Fig. 6A) was designed with an interactive panel (Fig. 6B) that contains pathway visualization checkboxes to select a type of enrichment method, a dropdown menu to select color code on the pathway, treatment checkboxes to compare the treatments. Pathview package version 1.30.1 (Luo & Brouwer, 2013) was used for the visualization of the pathways based on ORA or GSEA. An example of pathway is presented in Fig. 6C. SBGNview package (overlay omics data onto sbgn pathway diagrams) version 1.4.1 (Dong et al., 2021) was used to visualize the enriched pathway from ReactomePA (example is depicted in Fig. 6D). The plot_network and post_payload functions of STRINGdb package version 2.2.2 (Szklarczyk et al., 2019) were used to visualize the network of protein-protein interaction network (example is presented in Fig. 6E).

Literature information

A europepmc package version 0.4 was used to fetch abstracts, publication details, and bibliography metadata in the Literature info UI tab of the OMnalysis web application.

User manual

To help researchers, a user manual is drafted that provides step-by-step guidelines (Data S1).

Workflow for transcriptomics and proteomics data analysis

We performed ID conversion on transcriptomics and proteomics data taken as example matrices using Ensembl ID and uniport ID, respectively. OMnalysis is built to assign the most updated IDs. To the transcriptomic data matrix, the OMnalysis assigned 11,357 updated Ensembl IDs from a total of 11,398 uploaded Ensembl ID (Table S7). Next, to those updated IDs, OMnalysis was successful in assigning 10,951 human gene names, 10,932 HGNC symbols, 11,354 gene descriptions, and 7,463 UniProtKB/Swiss-prot ID (Table S7). In the case of proteomic data, 731 UniProt IDs were submitted to OMnalysis. Please note that, these 731 IDs include several repeated UniProt IDs from 4 treatments. From this list, 281 UniProt IDs were mapped to Ensembl gene ID, 277 to the gene name, 0 to HGNC symbol, 273 to gene description and 273 to UniProtKB/Swiss-Prot ID (Table S8).

Visualization, plotting and statistical filtering

When looking for the visualization, OMnalysis enables a user to generate separate plots as shown in Figs. 2–4 for each treatment. It generates PCA plots to identify the relationship and variability among the genes or proteins in the treatments. It produces scatter and volcano plots to visualize the significantly up and down regulated gene or proteins in each treatment. It makes Venn diagram to observe the intersection among genes or proteins and histogram to identify level of expression (logFC) of DEGs or DEPs in each treatment.

Further, OMnalysis provides the user functionality to filter out the genes or proteins which are not significantly expressed or identified. The software provides an option to set the cutoff values for logFC, logCPM, and P value for transcriptomic data (Fig. 4A), whereas logFC and P value for proteomic data. logFC ± 1.2, Log CPM >3, and P value <0.001 were set to select the DEGs from the transcriptomic data matrix. Whereas in proteomics logFC ±1.2 and P value <0.01 (FDR-adjusted P-value) were set to select the DEPs. The resultant list of the differentially expressed candidates was automatically transferred to the next level i.e., the functional and pathway analysis.

Functional and pathway analysis of example data

Gene ontology (GO) analysis were performed using ORA (at 0.05 q-value cutoff) and GSEA (0.5 P value cutoff) on the differentially expressed candidates selected above (Fig. S1). Using ORA (Table S9) and GSEA (Table S10), OMnalysis enriched the DEGs according to their biological process, molecular function, and cellular component. To minimize the chances of identification of false-positive GO IDs in ORA and GSEA, OMnalysis provides multiple hypothesis testing correction methods (Holm, Hochberg, Hommel, Bonferroni, Benjamini and Hochberg, BY and FDR) that enable adjustment of P value (Fig. S1). In OMnalysis, the user can select any single GO ID that is present in all the treatments and compare the expression of the gene with the heat maps. An example of the heatmap is presented in Fig. S3. In the same way, OMnalysis processed the set of the differentially expressed proteins (Fig. 4A and 4B) using ORA (0.05 q-value cutoff; Table S11) and GSEA (0.5 P value cutoff; Table S12) and enriched them according to the biological process, molecular function, and cellular component. OMnalysis can generate the same type of heatmaps as shown in Fig. S3 for proteomics.

For pathway analysis, OMnalysis used the list of differentially expressed candidates (DEGs and DEPs) and enriched them using ORA (Table S13) and GSEA (Table S14) against KEGG pathways (Fig. S2). For pathway enrichment analysis one of the multiple hypothesis testing correction methods can also be selected (Fig. S2). The network topology analysis (NTA) was also performed on OMnalysis to enrich the list of DEGs against various databases such as biocarta, panther, NCI, PharmGKB (Table S15). OMnalysis was also able to compare the list of differentially expressed candidates to the pathways available in Reactome (Table S16).

OMnalysis was extended to include the STRINGdb, which covers several databases (InterPro, SMART, PFAM domains, Reactome pathways, PubMed publications, UniProt Keywords, GO terms, and KEGG) to enrich the DEGs (Table S17).

One of the main features of OMnalysis is the visualization of the level of expression (using different color codes) of the given gene or protein in different treatments. As shown in Fig. S4 , the OMnalysis enriched DEGs from all four treatments to TNF signaling pathway (using ORA against KEGG pathways) and assigned color codes based on their expression values (logFC) in the gene box. This kind of representation gives a holistic view of the ongoing molecular activities simultaneously in all treatments in a given pathway. In the same way, OMnalysis used DEGs and visualized their level of elicitation on the enriched Reactome pathway (Fig. S5). To extend the functionality, OMnalysis visualizes the level of elicitation in protein-protein interaction network derived from STRINGdb in the form of halo color codes (Fig. S6).

In the case of the proteomics study, OMnalysis processed the list of DEPs using ORA (Table S18) and GSEA (Table S19) against the KEGG database. Whereas, the unavailability of supporting databases for bovine proteome hindered the pathway enrichment analysis using NTA and Reactome. While choosing STRINGdb, OMnalysis enriched the DEPs to InterPro, SMART, PFAM domains, Reactome pathways, PubMed publications, UniProt Keywords, GO terms, and KEGG (Table S20). For proteomics, OMnalysis can visualize the generated pathway in the same way as in Figs. S4–S6.

OMnalysis output formats

OMnalysis was developed considering the maximum flexibility and customization of the resulting output. Users were provided with the option to download the converted accession IDs, results from GO enrichment analysis and pathway enrichment analysis in CSV file format. All the diagrams produced in PCA, Plots, Statistical filtering, and GO heatmaps sections can be downloaded in TIFF, PNG, JPEG, and PDF format. In the case of the enriched pathway section, the visualization output can be downloaded in PNG image format. Furthermore, users have the option to adjust the dimension and resolution of the images generated in each tab of OMnalysis.