Example workflow
2023-04-25
Contents
This vignette shows an example workflow for ensemble biclustering analysis
with the mosbi package.
Every function of the package has a help page with a detailed documentation.
To access these type help(package=mosbi) in the R console.
0.0.1 Load packages
Import dependencies.
library(mosbi)0.0.2 Helper functions
Two additional functions are defined, to calculate z-scores of the data and to visualize the biclusters as a histogram.
z_score <- function(x, margin = 2) {
z_fun <- function(y) {
(y - mean(y, na.rm = TRUE)) / sd(y, na.rm = TRUE)
}
if (margin == 2) {
return(apply(x, margin, z_fun))
} else if (margin == 1) {
return(t(apply(x, margin, z_fun)))
}
}
bicluster_histo <- function(biclusters) {
cols <- mosbi::colhistogram(biclusters)
rows <- mosbi::rowhistogram(biclusters)
graphics::par(mfrow = c(1, 2))
hist(cols, main = "Column size ditribution")
hist(rows, main = "Row size ditribution")
}0.0.3 1. Download and prepare data
Biclustering will be done on a data matrix. As an example,
lipidomics dataset from the metabolights database will be used
https://www.ebi.ac.uk/metabolights/MTBLS562.
The data consists of 40 samples (columns) and 245 lipids (rows).
# get data
data(mouse_data)
mouse_data <- mouse_data[c(
grep(
"metabolite_identification",
colnames(mouse_data)
),
grep("^X", colnames(mouse_data))
)]
# Make data matrix
data_matrix <- z_score(log2(as.matrix(mouse_data[2:ncol(mouse_data)])), 1)
rownames(data_matrix) <- mouse_data$metabolite_identification
stats::heatmap(data_matrix)
The data has a gaussian-like distribution and no missing values,
so we can proceed with biclustering.
0.0.4 2. Compute biclusters
The mosbi package is able to work with results of different
biclustering algorithms. The approach unites the results
from different algorithms.
The results of four example algorithms will be computed and
converted to mosbi::bicluster objects. For a list of all
supported biclustering algorithms/packages type ?mosbi::get_biclusters
in the R console.
# Fabia
fb <- mosbi::run_fabia(data_matrix) # In case the algorithms throws an error,
#> Cycle: 0
Cycle: 20
Cycle: 40
Cycle: 60
Cycle: 80
Cycle: 100
Cycle: 120
Cycle: 140
Cycle: 160
Cycle: 180
Cycle: 200
Cycle: 220
Cycle: 240
Cycle: 260
Cycle: 280
Cycle: 300
Cycle: 320
Cycle: 340
Cycle: 360
Cycle: 380
Cycle: 400
Cycle: 420
Cycle: 440
Cycle: 460
Cycle: 480
Cycle: 500
# return an empty list
# isa2
BCisa <- mosbi::run_isa(data_matrix)
# Plaid
BCplaid <- mosbi::run_plaid(data_matrix)
#> layer: 0
#> 5882.744
#> layer: 1
#> [1] 0 107 10
#> [1] 1 101 10
#> [1] 30 101 10
#> [1] 31 39 10
#> [1] 32 39 9
#> [1] 33 37 9
#> [1] 34 37 9
#> [1] 35 37 9
#> [1] 60 37 9
#> [1] 2
#> [1] 154.413 0.000 0.000 0.000
#> back fitting 2 times
#> layer: 2
#> [1] 0 50 11
#> [1] 1 54 13
#> [1] 2 52 14
#> [1] 3 46 16
#> [1] 4 50 17
#> [1] 5 57 17
#> [1] 6 72 17
#> [1] 7 98 17
#> [1] 8 103 17
#> [1] 9 104 17
#> [1] 30 104 17
#> [1] 31 19 17
#> [1] 32 19 8
#> [1] 33 19 8
#> [1] 60 19 8
#> [1] 10
#> [1] 95.26758 0.00000 0.00000 0.00000
#> back fitting 2 times
#> layer: 3
#> [1] 0 111 13
#> [1] 1 95 13
#> [1] 2 90 13
#> [1] 3 84 13
#> [1] 4 81 13
#> [1] 5 79 13
#> [1] 6 78 13
#> [1] 30 78 13
#> [1] 31 5 13
#> [1] 32 5 5
#> [1] 33 5 5
#> [1] 60 5 5
#> [1] 7
#> [1] 16.69135 0.00000 0.00000 0.00000
#> back fitting 2 times
#> layer: 4
#> [1] 0 111 19
#> [1] 1 85 16
#> [1] 2 73 15
#> [1] 3 65 15
#> [1] 4 60 15
#> [1] 5 56 15
#> [1] 6 53 15
#> [1] 30 53 15
#> [1] 31 3 15
#> [1] 32 3 9
#> [1] 33 3 9
#> [1] 34 3 8
#> [1] 35 3 8
#> [1] 60 3 8
#> [1] 7
#> [1] 18.12064 0.00000 0.00000 0.00000
#> back fitting 2 times
#> layer: 5
#> [1] 0 40 6
#> [1] 1 63 7
#> [1] 2 81 7
#> [1] 3 89 7
#> [1] 4 89 6
#> [1] 5 94 6
#> [1] 30 94 6
#> [1] 31 33 6
#> [1] 32 33 6
#> [1] 33 31 6
#> [1] 34 31 6
#> [1] 35 30 6
#> [1] 36 30 6
#> [1] 37 30 6
#> [1] 60 30 6
#> [1] 6
#> [1] 212.2954 0.0000 0.0000 0.0000
#> back fitting 2 times
#> layer: 6
#> [1] 0 77 18
#> [1] 1 61 16
#> [1] 2 52 15
#> [1] 3 49 15
#> [1] 30 49 15
#> [1] 31 0 15
#> [1] 32
#> [1] 0 0 0 0
#>
#> Layer Rows Cols Df SS MS Convergence Rows Released Cols Released
#> 0 245 40 284 6239.96 21.97 NA NA NA
#> 1 37 9 45 273.49 6.08 1 64 1
#> 2 19 8 26 211.28 8.13 1 85 9
#> 3 5 5 9 27.52 3.06 1 73 8
#> 4 3 8 10 29.64 2.96 1 50 7
#> 5 30 6 35 327.88 9.37 1 64 0
# QUBIC
BCqubic <- mosbi::run_qubic(data_matrix)
# Merge results of all algorithms
all_bics <- c(fb, BCisa, BCplaid, BCqubic)
bicluster_histo(all_bics)
The histogram visualizes the distribution of bicluster sizes
(separately for the number of rows and columns of each bicluster).
The total number of found biclusters are given in the title.
0.0.5 3. Compute network
The next step of the ensemble approach is the computation of a
similarity network of biclusters. To filter for for similarities due to
random overlaps of biclusters, we apply an error
model (For more details refer to our publication).
Different similarity metrics are available.
For details type mosbi::bicluster_network in the R console.
bic_net <- mosbi::bicluster_network(all_bics, # List of biclusters
data_matrix, # Data matrix
n_randomizations = 5,
# Number of randomizations for the
# error model
MARGIN = "both",
# Use datapoints for metric evaluation
metric = 4, # Fowlkes–Mallows index
# For information about the metrics,
# visit the "Similarity metrics
# evaluation" vignette
n_steps = 1000,
# At how many steps should
# the cut-of is evaluated
plot_edge_dist = TRUE
# Plot the evaluation of cut-off estimation
)
#> Esimated cut-off: 0.05505506
The two resulting plot visualize the process of cut-off estimation.
The right plot show the remaining number of edges for the computed
bicluster network (red) and for randomizations of biclusters (black).
The vertical red line showed the threshold with the highest
signal-to-noise ratio (SNR). All evaluated SNRs are again
visualized in the left plot.
The next plot shows the bicluster similarity matrix. It reveals highly similar biclusters.
stats::heatmap(get_adjacency(bic_net))
0.0.5.1 Visualize network
Before the final step, extraction of bicluster communities (ensemble biclusters), the bicluster network can be layouted as a network.
plot(bic_net)
The networks are plotted using the igraph package. igraph
specific plotting parameters can be added.
For help type: ?igraph::plot.igraph
To see, which bicluster was generated by which algorithm, the following function can be executed:
mosbi::plot_algo_network(bic_net, all_bics, vertex.label = NA)
The downloaded data contains samples from different weeks of development. This can be visualized on the network, showing from which week the samples within a bicluster come from.
# Prepare groups for plotting
weeks <- vapply(
strsplit(colnames(data_matrix), "\\."),
function(x) {
return(x[1])
}, ""
)
names(weeks) <- colnames(data_matrix)
print(sort(unique(weeks))) # 5 colors required
#> [1] "X12W" "X24W" "X32W" "X4W" "X52W"
week_cols <- c("yellow", "orange", "red", "green", "brown")
# Plot network colored by week
mosbi::plot_piechart_bicluster_network(bic_net, all_bics, weeks,
week_cols,
vertex.label = NA
)
graphics::legend("topright",
legend = sort(unique(weeks)),
fill = week_cols, title = "Week"
)
Such a visualization is also possible for the samples:
# Prepare groups for plotting
samples <- vapply(
strsplit(colnames(data_matrix), "\\."),
function(x) {
return(x[2])
}, ""
)
names(samples) <- colnames(data_matrix)
samples_cols <- RColorBrewer::brewer.pal(
n = length(sort(unique(samples))),
name = "Set3"
)
# Plot network colored by week
mosbi::plot_piechart_bicluster_network(bic_net, all_bics, samples,
samples_cols,
vertex.label = NA
)
graphics::legend("topright",
legend = sort(unique(samples)),
fill = samples_cols, title = "Sample"
)
0.0.6 4. Extract louvain communities
Calculate the communities
coms <- mosbi::get_louvain_communities(bic_net,
min_size = 3,
bics = all_bics
)
# Only communities with a minimum size of 3 biclusters are saved.0.0.6.1 Visualization of the communities
# Plot all communities
for (i in seq(1, length(coms))) {
tmp_bics <- mosbi::select_biclusters_from_bicluster_network(
coms[[i]],
all_bics
)
mosbi::plot_piechart_bicluster_network(coms[[i]], tmp_bics,
weeks, week_cols,
main = paste0("Community ", i)
)
graphics::legend("topright",
legend = sort(unique(weeks)),
fill = week_cols, title = "Week"
)
cat("\nCommunity ", i, " conists of results from the
following algorithms:\n")
cat(get_bic_net_algorithms(coms[[i]]))
cat("\n")
}
#>
#> Community 1 conists of results from the
#> following algorithms:
#> fabia biclust-qubic
#>
#> Community 2 conists of results from the
#> following algorithms:
#> fabia isa2 biclust-qubic
#>
#> Community 3 conists of results from the
#> following algorithms:
#> isa2
#>
#> Community 4 conists of results from the
#> following algorithms:
#> isa2 biclust-qubic
#>
#> Community 5 conists of results from the
#> following algorithms:
#> isa2
#>
#> Community 6 conists of results from the
#> following algorithms:
#> isa2 biclust-plaid biclust-qubic0.0.6.2 Extraction of the communities
Finally, communities of the network can be extracted as ensemble biclusters.
The are saved as a list of mosbi::bicluster objects and therefore in the
same format as the imported results of all the algorithms.
With the parameters row_threshold & col_threshold,
the minimum occurrence of a row- or column-element in the biclusters
of a community can be defined.
ensemble_bicluster_list <- mosbi::ensemble_biclusters(coms, all_bics,
data_matrix,
row_threshold = .1,
col_threshold = .1
)1 Session Info
sessionInfo()
#> R version 4.3.0 RC (2023-04-13 r84269)
#> Platform: x86_64-pc-linux-gnu (64-bit)
#> Running under: Ubuntu 22.04.2 LTS
#>
#> Matrix products: default
#> BLAS: /home/biocbuild/bbs-3.17-bioc/R/lib/libRblas.so
#> LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.0
#>
#> locale:
#> [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
#> [3] LC_TIME=en_GB LC_COLLATE=C
#> [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
#> [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
#> [9] LC_ADDRESS=C LC_TELEPHONE=C
#> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
#>
#> time zone: America/New_York
#> tzcode source: system (glibc)
#>
#> attached base packages:
#> [1] stats graphics grDevices utils datasets methods base
#>
#> other attached packages:
#> [1] mosbi_1.6.0 BiocStyle_2.28.0
#>
#> loaded via a namespace (and not attached):
#> [1] tidyr_1.3.0 generics_0.1.3 sass_0.4.5
#> [4] utf8_1.2.3 class_7.3-21 lattice_0.21-8
#> [7] digest_0.6.31 magrittr_2.0.3 evaluate_0.20
#> [10] grid_4.3.0 RColorBrewer_1.1-3 bookdown_0.33
#> [13] fastmap_1.1.1 jsonlite_1.8.4 BiocManager_1.30.20
#> [16] purrr_1.0.1 fansi_1.0.4 scales_1.2.1
#> [19] modeltools_0.2-23 jquerylib_0.1.4 cli_3.6.1
#> [22] isa2_0.3.6 rlang_1.1.0 Biobase_2.60.0
#> [25] munsell_0.5.0 cachem_1.0.7 yaml_2.3.7
#> [28] tools_4.3.0 parallel_4.3.0 biclust_2.0.3
#> [31] dplyr_1.1.2 colorspace_2.1-0 ggplot2_3.4.2
#> [34] BiocGenerics_0.46.0 vctrs_0.6.2 R6_2.5.1
#> [37] stats4_4.3.0 lifecycle_1.0.3 magick_2.7.4
#> [40] QUBIC_1.28.0 MASS_7.3-59 pkgconfig_2.0.3
#> [43] RcppParallel_5.1.7 bslib_0.4.2 pillar_1.9.0
#> [46] gtable_0.3.3 glue_1.6.2 Rcpp_1.0.10
#> [49] tidyselect_1.2.0 xfun_0.39 tibble_3.2.1
#> [52] highr_0.10 flexclust_1.4-1 knitr_1.42
#> [55] fabia_2.46.0 igraph_1.4.2 htmltools_0.5.5
#> [58] rmarkdown_2.21 BH_1.81.0-1 compiler_4.3.0
#> [61] additivityTests_1.1-4.1