This vignette is made for users that are already familiar with the basic condiments workflow described in the first vignette. Here, we will show how to modify the default parameters for the first two steps of the workflow
# For analysis
library(condiments)
library(slingshot)
set.seed(21)We rely on the same toy dataset as the first vignette
data("toy_dataset", package = "condiments")
df <- toy_dataset$sd
rd <- as.matrix(df[, c("Dim1", "Dim2")])
sds <- slingshot(rd, df$cl)By default, the topologyTest function requires only
two inputs, the sds object and the condition
labels. To limit run time for the vignette, we also change the default
number of permutations used to generate trajectories under the null by
setting the rep argument to \(10\) instead of the default \(100\). As such, the test statistics might
be more variable.
top_res <- topologyTest(sds = sds, conditions = df$conditions, rep = 10)## Generating permuted trajectories## Running KS-mean testknitr::kable(top_res)| method | thresh | statistic | p.value |
|---|---|---|---|
| KS_mean | 0.01 | 0 | 1 |
The topologyTest function can be relatively slow on
large datasets. Moreover, when changing the method used to test the null
hypothesis that a common trajectory should be fitted, the first
permutation part of generating rep trajectories under the
null is identical. Therefore, we allow users to specify more than one
method and one value of the threshold. Here, we will use both the
Kolmogorov-Smirnov test test(Smirnov 1939)
and the classifier-test(Lopez-Paz and Oquab
2016).
top_res <- topologyTest(sds = sds, conditions = df$conditions, rep = 10,
methods = c("KS_mean", "Classifier"),
threshs = c(0, .01, .05, .1))## Generating permuted trajectories## Running KS-mean test## Running Classifier testknitr::kable(top_res)| method | thresh | statistic | p.value |
|---|---|---|---|
| KS_mean | 0 | 0.0070000 | 1.0000000 |
| KS_mean | 0.01 | 0.0000000 | 1.0000000 |
| KS_mean | 0.05 | 0.0000000 | 1.0000000 |
| KS_mean | 0.1 | 0.0000000 | 1.0000000 |
| Classifier | 0 | 0.4150000 | 0.9999821 |
| Classifier | 0.01 | 0.3800000 | 1.0000000 |
| Classifier | 0.05 | 0.3333333 | 1.0000000 |
| Classifier | 0.1 | 0.2833333 | 1.0000000 |
To see all methods avaible, use
/tmp/Rtmp8O3jlU/Rinst28def6a3ec799/condiments/help/topologyTest and look
at the methods argument.
For all methods but the KS test, additional paramters can be
specified, using a custom argument: args_classifier,
args_wass or args_mmd. See the help file for
given test more information on those parameters. For example, since the
default test based on the wasserstein distance and permutation test is
quite slow, we can pass a fast argument.
top_res <- topologyTest(sds = sds, conditions = df$conditions, rep = 10,
methods = "wasserstein_permutation",
args_wass = list(fast = TRUE, S = 100, iterations = 10^2))## Generating permuted trajectories## Running wassertsein permutation testknitr::kable(top_res)| method | thresh | statistic | p.value |
|---|---|---|---|
| wasserstein_permutation | NA | 1.356861 | 0.85 |
For now, the first part of the topologyTest has been designed for parallelisation using the BiocParallel package. For example, to run with 4 cores, you can run the following command
library(BiocParallel)
BPPARAM <- bpparam()
BPPARAM$progressbar <- TRUE
BPPARAM$workers <- 4
top_res <- topologyTest(sds = sds, conditions = df$conditions, rep = 100,
parallel = TRUE, BPPARAM = BPPARAM)
knitr::kable(top_res)The tests for the second test are much less compute-intensive, therefore there is no parallelisation. However, the other changes introduce in the previous section are still possible
prog_res <- progressionTest(sds, conditions = df$conditions)
knitr::kable(prog_res)| lineage | statistic | p.value |
|---|---|---|
| All | 5.504172 | 0 |
dif_res <- fateSelectionTest(sds, conditions = df$conditions)## note: only 1 unique complexity parameters in default grid. Truncating the grid to 1 .
##
## note: only 1 unique complexity parameters in default grid. Truncating the grid to 1 .knitr::kable(dif_res)| pair | statistic | p.value |
|---|---|---|
| 1vs2 | 0.6386486 | 1.9e-05 |
prog_res <- progressionTest(sds, conditions = df$conditions, method = "Classifier")
knitr::kable(prog_res)| lineage | statistic | p.value |
|---|---|---|
| All | 0.5890991 | 0.0043539 |
dif_res <- fateSelectionTest(sds, conditions = df$conditions, thresh = .05)## note: only 1 unique complexity parameters in default grid. Truncating the grid to 1 .
##
## note: only 1 unique complexity parameters in default grid. Truncating the grid to 1 .knitr::kable(dif_res)| pair | statistic | p.value |
|---|---|---|
| 1vs2 | 0.6121622 | 0.0004817 |
prog_res <- progressionTest(sds, conditions = df$conditions, method = "Classifier",
args_classifier = list(method = "rf"))## Warning in randomForest.default(x, y, mtry = param$mtry, ...): invalid mtry:
## reset to within valid range
## Warning in randomForest.default(x, y, mtry = param$mtry, ...): invalid mtry:
## reset to within valid range
## Warning in randomForest.default(x, y, mtry = param$mtry, ...): invalid mtry:
## reset to within valid range
## Warning in randomForest.default(x, y, mtry = param$mtry, ...): invalid mtry:
## reset to within valid range
## Warning in randomForest.default(x, y, mtry = param$mtry, ...): invalid mtry:
## reset to within valid range
## Warning in randomForest.default(x, y, mtry = param$mtry, ...): invalid mtry:
## reset to within valid range
## Warning in randomForest.default(x, y, mtry = param$mtry, ...): invalid mtry:
## reset to within valid range
## Warning in randomForest.default(x, y, mtry = param$mtry, ...): invalid mtry:
## reset to within valid range
## Warning in randomForest.default(x, y, mtry = param$mtry, ...): invalid mtry:
## reset to within valid range
## Warning in randomForest.default(x, y, mtry = param$mtry, ...): invalid mtry:
## reset to within valid range
## Warning in randomForest.default(x, y, mtry = param$mtry, ...): invalid mtry:
## reset to within valid range
## Warning in randomForest.default(x, y, mtry = param$mtry, ...): invalid mtry:
## reset to within valid range
## Warning in randomForest.default(x, y, mtry = param$mtry, ...): invalid mtry:
## reset to within valid range
## Warning in randomForest.default(x, y, mtry = param$mtry, ...): invalid mtry:
## reset to within valid range
## Warning in randomForest.default(x, y, mtry = param$mtry, ...): invalid mtry:
## reset to within valid range
## Warning in randomForest.default(x, y, mtry = param$mtry, ...): invalid mtry:
## reset to within valid range
## Warning in randomForest.default(x, y, mtry = param$mtry, ...): invalid mtry:
## reset to within valid range
## Warning in randomForest.default(x, y, mtry = param$mtry, ...): invalid mtry:
## reset to within valid range
## Warning in randomForest.default(x, y, mtry = param$mtry, ...): invalid mtry:
## reset to within valid range
## Warning in randomForest.default(x, y, mtry = param$mtry, ...): invalid mtry:
## reset to within valid range
## Warning in randomForest.default(x, y, mtry = param$mtry, ...): invalid mtry:
## reset to within valid range
## Warning in randomForest.default(x, y, mtry = param$mtry, ...): invalid mtry:
## reset to within valid range## note: only 1 unique complexity parameters in default grid. Truncating the grid to 1 .knitr::kable(prog_res)| lineage | statistic | p.value |
|---|---|---|
| All | 0.517027 | 0.3192952 |
dif_res <- fateSelectionTest(sds, conditions = df$conditions)## note: only 1 unique complexity parameters in default grid. Truncating the grid to 1 .
##
## note: only 1 unique complexity parameters in default grid. Truncating the grid to 1 .knitr::kable(dif_res)| pair | statistic | p.value |
|---|---|---|
| 1vs2 | 0.6611712 | 8e-07 |
For all of the above procedures, it is important to note that we are making multiple comparisons. The p-values we obtain from these tests should be corrected for multiple testing, especially for trajectories with a large number of lineages.
That said, trajectory inference is often one of the last computational methods in a very long analysis pipeline (generally including gene-level quantification, gene filtering / feature selection, and dimensionality reduction). Hence, we strongly discourage the reader from putting too much faith in any p-value that comes out of this analysis. Such values may be useful suggestions, indicating particular features or cells for follow-up study, but should generally not be treated as meaningful statistical quantities.
If some commands and parameters are still unclear after going through this vignette, do not hesitate to open an issue on the condiments Github repository.
sessionInfo()## R Under development (unstable) (2024-10-21 r87258)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 24.04.1 LTS
##
## Matrix products: default
## BLAS: /home/biocbuild/bbs-3.21-bioc/R/lib/libRblas.so
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.12.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] stats4 stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] caret_6.0-94 lattice_0.22-6
## [3] viridis_0.6.5 viridisLite_0.4.2
## [5] RColorBrewer_1.1-3 ggplot2_3.5.1
## [7] tidyr_1.3.1 dplyr_1.1.4
## [9] slingshot_2.15.0 TrajectoryUtils_1.15.0
## [11] SingleCellExperiment_1.29.1 SummarizedExperiment_1.37.0
## [13] Biobase_2.67.0 GenomicRanges_1.59.1
## [15] GenomeInfoDb_1.43.2 IRanges_2.41.2
## [17] S4Vectors_0.45.2 BiocGenerics_0.53.3
## [19] generics_0.1.3 MatrixGenerics_1.19.0
## [21] matrixStats_1.4.1 princurve_2.1.6
## [23] condiments_1.15.1 knitr_1.49
##
## loaded via a namespace (and not attached):
## [1] jsonlite_1.8.9 magrittr_2.0.3
## [3] spatstat.utils_3.1-1 ggbeeswarm_0.7.2
## [5] farver_2.1.2 rmarkdown_2.29
## [7] zlibbioc_1.53.0 vctrs_0.6.5
## [9] DelayedMatrixStats_1.29.0 htmltools_0.5.8.1
## [11] S4Arrays_1.7.1 BiocNeighbors_2.1.1
## [13] SparseArray_1.7.2 pROC_1.18.5
## [15] sass_0.4.9 parallelly_1.40.1
## [17] bslib_0.8.0 plyr_1.8.9
## [19] lubridate_1.9.3 cachem_1.1.0
## [21] igraph_2.1.1 lifecycle_1.0.4
## [23] iterators_1.0.14 pkgconfig_2.0.3
## [25] rsvd_1.0.5 Matrix_1.7-1
## [27] R6_2.5.1 fastmap_1.2.0
## [29] GenomeInfoDbData_1.2.13 future_1.34.0
## [31] digest_0.6.37 colorspace_2.1-1
## [33] scater_1.35.0 irlba_2.3.5.1
## [35] beachmat_2.23.3 labeling_0.4.3
## [37] randomForest_4.7-1.2 fansi_1.0.6
## [39] timechange_0.3.0 mgcv_1.9-1
## [41] httr_1.4.7 abind_1.4-8
## [43] compiler_4.5.0 rngtools_1.5.2
## [45] proxy_0.4-27 withr_3.0.2
## [47] doParallel_1.0.17 BiocParallel_1.41.0
## [49] MASS_7.3-61 lava_1.8.0
## [51] DelayedArray_0.33.3 ModelMetrics_1.2.2.2
## [53] tools_4.5.0 vipor_0.4.7
## [55] beeswarm_0.4.0 future.apply_1.11.3
## [57] nnet_7.3-19 glue_1.8.0
## [59] nlme_3.1-166 grid_4.5.0
## [61] reshape2_1.4.4 recipes_1.1.0
## [63] gtable_0.3.6 class_7.3-22
## [65] data.table_1.16.2 BiocSingular_1.23.0
## [67] ScaledMatrix_1.15.0 utf8_1.2.4
## [69] XVector_0.47.0 ggrepel_0.9.6
## [71] RANN_2.6.2 foreach_1.5.2
## [73] pillar_1.9.0 stringr_1.5.1
## [75] limma_3.63.2 Ecume_0.9.2
## [77] splines_4.5.0 survival_3.7-0
## [79] tidyselect_1.2.1 scuttle_1.17.0
## [81] pbapply_1.7-2 transport_0.15-4
## [83] gridExtra_2.3 xfun_0.49
## [85] statmod_1.5.0 hardhat_1.4.0
## [87] distinct_1.19.0 timeDate_4041.110
## [89] stringi_1.8.4 UCSC.utils_1.3.0
## [91] yaml_2.3.10 evaluate_1.0.1
## [93] codetools_0.2-20 kernlab_0.9-33
## [95] tibble_3.2.1 cli_3.6.3
## [97] rpart_4.1.23 munsell_0.5.1
## [99] jquerylib_0.1.4 Rcpp_1.0.13-1
## [101] globals_0.16.3 spatstat.univar_3.1-1
## [103] parallel_4.5.0 gower_1.0.1
## [105] doRNG_1.8.6 sparseMatrixStats_1.19.0
## [107] listenv_0.9.1 ipred_0.9-15
## [109] scales_1.3.0 prodlim_2024.06.25
## [111] e1071_1.7-16 purrr_1.0.2
## [113] crayon_1.5.3 rlang_1.1.4