bulkAnalyseR: An accessible, interactive pipeline for analysing and sharing bulk sequencing results

Loading an expression matrix

For this vignette we are using a subset of the count matrix for an experiment included in a 2019 paper by Yang et al. Rows represent genes/features and columns represent samples:

counts.in <- system.file("extdata", "expression_matrix.csv", package = "bulkAnalyseR")
exp <- as.matrix(read.csv(counts.in, row.names = 1))
head(exp)
##>                    SRR7624365 SRR7624366 SRR7624371 SRR7624372 SRR7624375
##> ENSMUSG00000102693          2          0          0          0          0
##> ENSMUSG00000051951          6          4          2          0         47
##> ENSMUSG00000102851          0          0          0          0          0
##> ENSMUSG00000103377          0          0          2          0          8
##> ENSMUSG00000104017          0          0          0          2          0
##> ENSMUSG00000103025          0          0          0          0          0
##>                    SRR7624376
##> ENSMUSG00000102693          0
##> ENSMUSG00000051951         37
##> ENSMUSG00000102851          1
##> ENSMUSG00000103377          6
##> ENSMUSG00000104017          2
##> ENSMUSG00000103025          3

Defining metadata

The subset of the dataset we are using has samples from 3 timepoints: 0h, 12h and 36h, each with 2 biological replicates. We define a metadata table detailing which sample correspond to which timepoint:

meta <- data.frame(
  srr = colnames(exp), 
  timepoint = rep(c("0h", "12h", "36h"), each = 2)
)
meta
##>          srr timepoint
##> 1 SRR7624365        0h
##> 2 SRR7624366        0h
##> 3 SRR7624371       12h
##> 4 SRR7624372       12h
##> 5 SRR7624375       36h
##> 6 SRR7624376       36h

This metadata table should be a data frame containing at minimum two columns: the first column must contain the column names of the expression.matrix, while the last column is assumed to contain the experimental conditions that will be tested for differential expression.

Denoising and normalisation

Before using the expression matrix to create our shiny app, some preprocessing should be performed. bulkAnalyseR contains the function preprocessExpressionMatrix which takes the expression matrix as input then denoises the data using noisyR and normalises it (quantile normalisation by default, other methods can bespecified using the normalisation.method parameter). By specifying output.plot = TRUE, you can also print the expression-similarity line plots from noisyR to console and you can specify further parameters from the noisyR noisyr_counts.

exp.proc <- preprocessExpressionMatrix(exp, output.plot = TRUE)
##> >>> noisyR counts approach pipeline <<<
##> The input matrix has 38705 rows and 6 cols
##>     number of genes: 38705
##>     number of samples: 6
##> Calculating the number of elements per window
##>     the number of elements per window is 3870
##>     the step size is 193
##>     the selected similarity metric is correlation_pearson
##>   Working with sample 1
##>   Working with sample 2
##>   Working with sample 3
##>   Working with sample 4
##>   Working with sample 5
##>   Working with sample 6
##> Calculating noise thresholds for 6 samples...
##>     similarity.threshold = 0.25
##>     method.chosen = Boxplot-IQR
##> Denoising expression matrix...
##>     removing noisy genes
##>     adjusting matrix
##> >>> Done! <<<
##> Performing quantile normalisation...
##> Done!

It is not recommended to use data which has not been denoised and normalised as input to generateShinyApp. You can also perform your own preprocessing outside preprocessExpressionMatrix.

Home

This tab contains some useful links (including this vignette!) to get you started with bulkAnalyseR.

Sample Select

This tab allows you to select a subset of samples to use for the analysis. In particular this allows you to exclude low quality samples (which you could discover using QC tab) or focus on particular parts of the data. Once you have selected the samples you want, make sure to press the ‘Use the selected sample!’ button to confirm your choice. If you do not use this button or select samples then all samples will be used.

Quality check (QC)

The Quality Check (QC) tab includes a Jaccard Similarity Index (JSI) heatmap, and a PCA dimensionality reduction, with groups based on the metadata information. This enables a high-level overview of the similarity across samples, usually reflecting the experimental design quite closely and serving as a sanity check. Within this panel, there is also the option to show the MA plot between any 2 columns - this can be used to further investigate any sample similarity or differences found through the JSI and PCA plots.

jaccard_heatmap(
  expression.matrix = exp.proc,
  metadata = meta,
  top.annotation.ids = 2,
  n.abundant = 500, 
  show.values = FALSE,
  show.row.column.names = (nrow(meta) <= 20)
)

plot_pca(
  expression.matrix = exp.proc,
  metadata = meta,
  annotation.id = 2,
  n.abundant = 500,
  show.labels = TRUE,
  show.ellipses = TRUE
)

Differential expression (DE)

The differential expression (DE) tab can be used to compare any of the conditions from the columns of the metadata table using the edgeR or DESeq2 packages with customisable parameters (log2 fold change and adjusted p-value threshold). The output is a list of differentially expressed genes along with log2 fold changes and adjusted p-values. Within the app, the differentially expressed genes can be explored using an interactive table (sorted by absolute log2 fold change) and this can also be saved as a comma-separated file.

# first create annotation table
anno <- AnnotationDbi::select(
  org.Mm.eg.db::org.Mm.eg.db,
  keys = rownames(exp.proc),
  keytype = 'ENSEMBL',
  columns = 'SYMBOL'
) %>%
  distinct(ENSEMBL, .keep_all = TRUE) %>%
  mutate(NAME = ifelse(is.na(SYMBOL), ENSEMBL, SYMBOL))

DEtable <- DEanalysis_edger(
  expression.matrix = exp.proc[, 1:4],
  condition = meta$timepoint[1:4],
  var1 = "0h",
  var2 = "12h",
  anno = anno
)
DEtable_deseq <- DEanalysis_deseq2(
  expression.matrix = exp.proc[, 1:4],
  condition = meta$timepoint[1:4],
  var1 = "0h",
  var2 = "12h",
  anno = anno
)

In this tab, you can also select genes of interest by clicking on the table, which can be used in the DE summary and Volcano/MA tabs. There is also a button to choose the top 50 DE genes (by absolute log2 fold change) or reset the selection.

Volcano/MA plots

Assuming differential expression has been performed, the Volcano/MA panel then visualises the DE genes, allowing interactive exploration of the data by searching for specific genes of interest or using the selected top DE genes and clicking on genes on the plot itself to gain more information and generate hypotheses.

volcano_plot(
  genes.de.results = DEtable,
  pval.threshold = 0.05, 
  lfc.threshold = 1,
  log10pval.cap = TRUE
)

ma_plot(
  genes.de.results = DEtable,
  pval.threshold = 0.05, 
  lfc.threshold = 1
)

DE summary

The DE summary panel allows further visualisation using the genes called DE or selected by the user in the DE panel.

genes_heatmap <- filter(arrange(DEtable, desc(abs(log2FC))), pvalAdj < 0.05)$gene_id
expression_heatmap(
  expression.matrix.subset = exp.proc[genes_heatmap, ],
  top.annotation.ids = 2,
  metadata = meta,
  type = "Z-score",
  show.column.names = (nrow(meta) <= 20)
)

genes_de <- filter(DEtable, pvalAdj < 0.05, abs(log2FC) > 1)$gene_id
plot_pca(
  expression.matrix = exp.proc[genes_de, ],
  metadata = meta,
  annotation.id = 2,
  n.abundant = NULL,
  show.labels = TRUE,
  show.ellipses = TRUE
)

Enrichment

The enrichment tab uses the differentially expressed genes to run an enrichment analysis using g:profiler with a user selected output of the standard GO terms, the KEGG and reactome pathway databases, and miRNA and transcription factor lists. The plot shows the log10 BH-corrected p-value for each term, separated by category. The background set used is all the genes in the expression matrix - if noisyR preprocessing is used then this will be all the genes above the s/n threshold.

The output is both the results table from g:profiler which can be downloaded as well as a jitter plot which can be clicked to identify the most significant terms.

gostres <- gprofiler2::gost(
  query = genes_de,
  organism = "mmusculus",
  correction_method = 'fdr',
  custom_bg = DEtable$gene_id,
  evcodes = TRUE
)
gostres$result <-  mutate(
  gostres$result,
  parents = sapply(.data$parents, toString),
  intersection_names = sapply(.data$intersection, function(x){
    ensids <- strsplit(x, split = ",")[[1]]
    names <- DEtable$gene_name[match(ensids, DEtable$gene_id)]
    paste(names, collapse = ",")
  })
)

Expression patterns

The expression pattern tab allows the creation of expression patterns to identify potential genes of interest across a variety of conditions. The most common application of this is a time series, but it could be suitable for another logical progression between conditions. To define the series, the user must select a column of the metadata and drag states into the “Series of states to use” area.

The pattern identification is done by calculating a confidence interval for each gene in each condition, using all samples in that condition and the number of standard deviations away from the mean provided. The pattern between two consecutive conditions is defined as straight (S) if the intervals overlap and up (U) or down (D) if they don’t. The full expression pattern is the concatenation of individual patterns (for example, “UUS” for 4 conditions).

The grouped expression matrix can then be downloaded, showing which pattern each gene was assigned to. Plots are also created for the genes in the selected pattern (“Pattern to plot”).

condition <- factor(meta$timepoint)
tbl <- calculate_condition_mean_sd_per_gene(exp.proc, condition)
patterns <- make_pattern_matrix(tbl, n_sd = 2)[, "pattern"]

Cross plots

The cross plot tab allows you to compare two differential expression analyses against each other, for example two comparisons of interest or the same comparison using edgeR and DESeq2. The plot shows the log2 fold change of the two differential expression calls on each axis. Genes which are DE in both comparisons are coloured purple, in comparison 1 but not comparison 2 in blue and in comparison 2 but not comparison 1 in red. You can label selected genes and click on genes on the plot itself to gain more information and generate hypotheses.

cross_plot(
  DEtable1 = DEtable,
  DEtable2 = DEtable_deseq,
  DEtable1Subset = filter(DEtable, pvalAdj < 0.05, abs(log2FC) > 1),
  DEtable2Subset = filter(DEtable_deseq, pvalAdj < 0.05, abs(log2FC) > 1),
  lfc.threshold = 1
)

Gene regulatory networks (GRN)

bulkAnalyseR enables the creation of small gene regulatory networks (GRNs) to facilitate further exploration and hypothesis generation based on genes of interest. Target genes can be selected and a small network with them as targets can be generated by clicking the “Start GRN inference” button. Within a bulkAnalyseR app, there are 2 separate tabs where inferred GRNs can be calculated and examined. Firstly, in GRN inference, up to four different GRNs can be inferred at once, to investigate differences between experimental conditions. If more than one GRN is generated, regulators in common are coloured differently and an UpSet plot is shown, summarising the commonalities of edges between the networks. Secondly, in the DE & enrichment GRN tab, a GRN can be inferred using only differentially expressed genes (using the comparison selected in the DE tab) and additional information relating these genes to transcription factors and/or miRNAs using results from the enrichment tab can be appended to the resulting network. In both cases, the number of regulators can then be adjusted and the network plots downloaded in interactive html format. Note that target genes can also regulate each other if the selected genes are functionally similar.

infer_GRN(
  expression.matrix = exp.proc, 
  metadata = meta, 
  anno = anno,
  seed = 13,
  targets = "Trpm1", 
  condition = "timepoint", 
  samples = c("0h", "12h"), 
  inference_method = "GENIE3"
)

Adding extra panels

Alongside the default panels, you can also define your own panels and add them to the app. As an example, we could add an extra QC panel (in this case it will be exactly the same). Note that extra panels appear at the top level, and all parameters are converted to lists by generateShinyApp and thus need to be supplied as such. It might be useful to explore the app.R file generated by the command below to understand the structure:

generateShinyApp(
  shiny.dir = "shiny_Yang2019_ExtraQC",
  app.title = "Shiny app for visualisation of three timepoints from the Yang 2019 data - extra QC",
  modality = "RNA",
  expression.matrix = exp.proc,
  metadata = meta,
  organism = "mmusculus",
  org.db = "org.Mm.eg.db",
  panels.extra = tibble::tibble(
    name = "RNA2",
    UIfun = "modalityPanelUI", 
    UIvars = "'RNA2', metadata[[1]], NA, 'QC'", 
    serverFun = "modalityPanelServer", 
    serverVars = "'RNA2', expression.matrix[[1]], metadata[[1]], anno[[1]], NA, 'QC'"
  )
)

If you need to add extra data or package imports for the extra panel(s) then you can do this using the data.extra and packages.extra parameters. Make sure you have the extra data loaded when you create the app. Note that the data.exta object will be available in the app as a single list. For example:

extra.data1 = matrix(rnorm(36),nrow=6)
extra.data2 = matrix(rnorm(60),nrow=10)

generateShinyApp(
  shiny.dir = "shiny_Yang2019_ExtraQC",
  app.title = "Shiny app for visualisation of three timepoints from the Yang 2019 data - extra QC",
  modality = "RNA",
  expression.matrix = exp.proc,
  metadata = meta,
  organism = "mmusculus",
  org.db = "org.Mm.eg.db",
  panels.extra = tibble::tibble(
    name = "RNA2",
    UIfun = "modalityPanelUI", 
    UIvars = "'RNA2', metadata[[1]], NA, 'QC'", 
    serverFun = "modalityPanelServer", 
    serverVars = "'RNA2', expression.matrix[[1]], metadata[[1]], anno[[1]], NA, 'QC'"
  ),
  data.extra = list(extra.data1, extra.data2),
  packages.extra = "somePackage",
)

sessionInfo()
##> R version 4.2.2 (2022-10-31 ucrt)
##> Platform: x86_64-w64-mingw32/x64 (64-bit)
##> Running under: Windows 10 x64 (build 22621)
##> 
##> Matrix products: default
##> 
##> locale:
##> [1] LC_COLLATE=C                           
##> [2] LC_CTYPE=English_United Kingdom.utf8   
##> [3] LC_MONETARY=English_United Kingdom.utf8
##> [4] LC_NUMERIC=C                           
##> [5] LC_TIME=English_United Kingdom.utf8    
##> 
##> attached base packages:
##> [1] stats     graphics  grDevices utils     datasets  methods   base     
##> 
##> other attached packages:
##> [1] ggplot2_3.4.0      dplyr_1.0.10       bulkAnalyseR_1.1.0
##> 
##> loaded via a namespace (and not attached):
##>  [1] tidyselect_1.2.0      xfun_0.35             bslib_0.4.1          
##>  [4] lattice_0.20-45       splines_4.2.2         colorspace_2.0-3     
##>  [7] vctrs_0.5.1           generics_0.1.3        htmltools_0.5.4      
##> [10] yaml_2.3.6            mgcv_1.8-41           utf8_1.2.2           
##> [13] noisyr_1.0.0          rlang_1.0.6           jquerylib_0.1.4      
##> [16] pillar_1.8.1          later_1.3.0           glue_1.6.2           
##> [19] withr_2.5.0           DBI_1.1.3             foreach_1.5.2        
##> [22] lifecycle_1.0.3       stringr_1.5.0         munsell_0.5.0        
##> [25] gtable_0.3.1          codetools_0.2-18      evaluate_0.19        
##> [28] labeling_0.4.2        knitr_1.41            fastmap_1.1.0        
##> [31] httpuv_1.6.7          fansi_1.0.3           highr_0.9            
##> [34] preprocessCore_1.60.0 Rcpp_1.0.9            xtable_1.8-4         
##> [37] scales_1.2.1          promises_1.2.0.1      cachem_1.0.6         
##> [40] jsonlite_1.8.4        farver_2.1.1          mime_0.12            
##> [43] digest_0.6.31         stringi_1.7.8         shiny_1.7.3          
##> [46] grid_4.2.2            cli_3.4.1             tools_4.2.2          
##> [49] magrittr_2.0.3        philentropy_0.7.0     sass_0.4.4           
##> [52] tibble_3.1.8          pkgconfig_2.0.3       Matrix_1.5-1         
##> [55] ellipsis_0.3.2        assertthat_0.2.1      rmarkdown_2.18       
##> [58] rstudioapi_0.14       iterators_1.0.14      R6_2.5.1             
##> [61] nlme_3.1-160          compiler_4.2.2