VSClust with Bioconductor objects
Introduction
For a more detailed explanation of the VSClust function and the workflow, please take a look on the vignette for running the VSClust workflow.
Here, we present an example script to integrate the clustering with
data object from Bioconductor, such as QFeatures,
SummarizedExperiment and
MultiAssayExperiment.
Installation and additional packages
Use the common Bioconductor commands for installation:
# uncomment in case you have not installed vsclust yet
#if (!require("BiocManager", quietly = TRUE))
# install.packages("BiocManager")
#BiocManager::install("vsclust")The full functionality can be obtained by additionally installing and
loading the packages yaml, shiny,
clusterProfiler, and matrixStats.
Initialization
Here, we define the different parameters for the data set
RNASeq2GeneNorm from the miniACC object.
The number of replicates and experimental conditions will be retrieved automatically by specifying the metadata for the grouping.
#### Input parameters, only read when now parameter file was provided #####
## All principal parameters for running VSClust can be defined as in the
## shiny app at computproteomics.bmb.sdu.dk/Apps/VSClust
# name of study
Experiment <- "miniACC"
# Paired or unpaired statistical tests when carrying out LIMMA for
# statistical testing
isPaired <- FALSE
# Number of threads to accelerate the calculation (use 1 in doubt)
cores <- 1
# If 0 (default), then automatically estimate the cluster number for the
# vsclust run from the Minimum Centroid Distance
PreSetNumClustVSClust <- 0
# If 0 (default), then automatically estimate the cluster number for the
# original fuzzy c-means from the Minimum Centroid Distance
PreSetNumClustStand <- 0
# max. number of clusters when estimating the number of clusters.
# Higher numbers can drastically extend the computation time.
maxClust <- 10 Statistics and data preprocessing
At first, we load will log-transform the original data and normalize it to the median. Statistical testing will be applied on the resulting object. After estimating the standard deviations, the matrix consists of the averaged quantitative feature values and a last column for the standard deviations of the features.
We will separate the samples according to their
OncoSign.
data(miniACC, package="MultiAssayExperiment")
# log-transformation and remove of -Inf values
logminiACC <- log2(assays(miniACC)$RNASeq2GeneNorm)
logminiACC[!is.finite(logminiACC)] <- NA
# normalize to median
logminiACC <- t(t(logminiACC) - apply(logminiACC, 2, median, na.rm=TRUE))
miniACC2 <- c(miniACC, log2rnaseq = logminiACC, mapFrom=1L)## Warning: Assuming column order in the data provided
## matches the order in 'mapFrom' experiment(s) colnames
#### running statistical analysis and estimation of individual variances
statOut <- PrepareSEForVSClust(miniACC2, "log2rnaseq",
coldatname = "OncoSign",
isPaired=isPaired, isStat=TRUE)## -- The following categories will be used as experimental
## conditions:
## CN2
## TP53/NF1
## TERT/ZNRF3
## CN1
## Unclassified
## NA
## CTNNB1## -- Extracted NumReps: 21 and NumCond: 7
We can see that there is no good separation of cancer signatures on the PCA plot.
Estimation of cluster number
There is no simple way to find the optimal number of clusters in a data set. For obtaining this number, we run the clustering for different cluster numbers and evaluate them via so-called validity indices, which provide information about suitable cluster numbers. VSClust uses mainly the “Maximum centroid distances” that denotes the shortest distance between any of the centroids. Alternatively, one can inspect the Xie Beni index.
The output of estimClustNum contains the suggestion for
the number of clusters.
We further visualize the outcome.
#### Estimate number of clusters with maxClust as maximum number clusters to run
#### the estimation with
ClustInd <- estimClustNum(statOut$dat, maxClust=maxClust, cores=cores)## Running cluster number 3## Running cluster number 4## Running cluster number 5## Running cluster number 6## Running cluster number 7## Running cluster number 8## Running cluster number 9## Running cluster number 10#### Use estimate cluster number or use own
if (PreSetNumClustVSClust == 0)
PreSetNumClustVSClust <- optimalClustNum(ClustInd)
if (PreSetNumClustStand == 0)
PreSetNumClustStand <- optimalClustNum(ClustInd, method="FCM")
#### Visualize
estimClust.plot(ClustInd)
Both validity indices agree with each other and suggest 7 as the most
reasonable estimate for the cluster number. However, we can also see
that this decreases the number of clustered features quite drastically
from over 150 to about 90.
Run final clustering
Now we run the clustering again with the optimal parameters from the estimation. One can take alternative numbers of clusters corresponding to large decays in the Minimum Centroid Distance or low values of the Xie Beni index.
First, we carry out the variance-sensitive method
#### Run clustering (VSClust and standard fcm clustering
ClustOut <- runClustWrapper(statOut$dat,
PreSetNumClustVSClust, NULL,
VSClust=TRUE,
cores=cores)
Bestcl <- ClustOut$Bestcl
VSClust_cl <- Bestcl
We see how different groups of genes show distinctive pattern of their expression in different oncological signatures.
## R version 4.4.2 (2024-10-31)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 24.04.1 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.26.so; LAPACK version 3.12.0
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 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: Etc/UTC
## tzcode source: system (glibc)
##
## attached base packages:
## [1] stats4 stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] MultiAssayExperiment_1.33.0 SummarizedExperiment_1.37.0
## [3] Biobase_2.67.0 GenomicRanges_1.59.0
## [5] GenomeInfoDb_1.43.0 IRanges_2.41.0
## [7] S4Vectors_0.45.1 BiocGenerics_0.53.2
## [9] generics_0.1.3 MatrixGenerics_1.19.0
## [11] matrixStats_1.4.1 vsclust_1.9.3
## [13] BiocStyle_2.35.0
##
## loaded via a namespace (and not attached):
## [1] gtable_0.3.6 xfun_0.49 bslib_0.8.0
## [4] ggplot2_3.5.1 lattice_0.22-6 vctrs_0.6.5
## [7] tools_4.4.2 parallel_4.4.2 tibble_3.2.1
## [10] fansi_1.0.6 BiocBaseUtils_1.9.0 pkgconfig_2.0.3
## [13] Matrix_1.7-1 lifecycle_1.0.4 GenomeInfoDbData_1.2.13
## [16] stringr_1.5.1 compiler_4.4.2 statmod_1.5.0
## [19] munsell_0.5.1 httpuv_1.6.15 htmltools_0.5.8.1
## [22] sys_3.4.3 buildtools_1.0.0 sass_0.4.9
## [25] yaml_2.3.10 later_1.3.2 pillar_1.9.0
## [28] crayon_1.5.3 jquerylib_0.1.4 DelayedArray_0.33.1
## [31] cachem_1.1.0 limma_3.63.2 abind_1.4-8
## [34] mime_0.12 tidyselect_1.2.1 digest_0.6.37
## [37] stringi_1.8.4 reshape2_1.4.4 dplyr_1.1.4
## [40] splines_4.4.2 maketools_1.3.1 fastmap_1.2.0
## [43] grid_4.4.2 colorspace_2.1-1 cli_3.6.3
## [46] SparseArray_1.7.1 magrittr_2.0.3 S4Arrays_1.7.1
## [49] utf8_1.2.4 promises_1.3.0 UCSC.utils_1.3.0
## [52] scales_1.3.0 rmarkdown_2.29 XVector_0.47.0
## [55] httr_1.4.7 qvalue_2.39.0 shiny_1.9.1
## [58] evaluate_1.0.1 knitr_1.49 rlang_1.1.4
## [61] Rcpp_1.0.13-1 xtable_1.8-4 glue_1.8.0
## [64] BiocManager_1.30.25 jsonlite_1.8.9 plyr_1.8.9
## [67] R6_2.5.1 zlibbioc_1.52.0