# Load required packages
suppressPackageStartupMessages({
library(treekoR)
library(SingleCellExperiment)
library(ggtree)
})# Install the development version from GitHub:
# install.packages("devtools")
devtools::install_github("adam2o1o/treekoR")
library(treekoR)treekoR is a novel framework that aims to utilise the hierarchical nature of single cell cytometry data, to find robust and interpretable associations between cell subsets and patient clinical end points. This is achieved by deriving the tree structure of cell clusters, followed by measuring the %parent (proportions of each node in the tree relative to the number of cells belonging to the immediate parent node), in addition to the %total (proportion of cells in each node relative to all cells). These proportions are then used in significance testing and classification models to determine which cell subpopulation proportions most correlated with the patient clinical outcome of interest. treekoR then provides an interactive visualisation which helps to highlight these results.
SingleCellExperiment containing samples of flow cytometry expression data from 39 patients. This data represents a subset of a dataset that was originally used by De Biasi et al. (2020) for the characterisation of CD8+ T cells, comparing between COVID-19 patients and healthy controls.data(COVIDSampleData)
sce <- DeBiasi_COVID_CD8_samptreekoR requires the following information in the variables:
exprs: Single cell expression data (\(n \times p\)), where \(p\) is the number of markers, and \(n\) is the number of cellsclusters: a vector of length \(n\) representing the cell type or cluster of each cell (can be character or numeric)classes: a vector of length \(n\) containing the patient outcome/class each cell belongs tosamples: a vector of length \(n\) identifying the patient each cell belongs toIn this example: the clusters contain 100 clusters generated by FlowSOM; classes identify whether the cell belongs to a COVID-19 or healthy patient; and samples identify which cell the patient comes from.
exprs <- t(assay(sce, "exprs"))
clusters <- colData(sce)$cluster_id
classes <- colData(sce)$condition
samples <- colData(sce)$sample_idThe scaled median marker expression for each cluster is calculated which is used to construct a hierarchical tree.
In this step, the choice of hierarchical aggregation method (which determines the structure of the tree) is determined. By default the framework chooses HOPACH to construct the tree via the hierarchy_method argument, however any of the methods in hclust can be used (see 3.5.1).
clust_tree <- getClusterTree(exprs,
clusters,
hierarchy_method="hopach")Proportions of each cell cluster in the tree are calculated - both the proportion relative to all and proportion relative to the hierarchical parent. These proportions are used in a two sample t-test, testing for equal means between the patient clinical outcome using both types of proportions.
tested_tree <- testTree(phylo=clust_tree$clust_tree,
clusters=clusters,
samples=samples,
classes=classes,
pos_class_name=NULL)node: unique identifier for each node in the hierarchical treeparent: the node of the parentisTip: whether the node is a leaf node in the treeclusters: the clusters belonging to the corresponding nodestat_all: test statistic obtained from testing between conditions using the proportion of the node relative to all cells (%total) in each sample. pval_total is the corresponding p-value (unadjusted)stat_parent: test statistic obtained from testing between conditions using the proportion of the node relative to cells in the parent node (%parent) in each sample. pval_parent is the corresponding p-value (unadjusted)res_df <- getTreeResults(tested_tree)
head(res_df, 10)
#> parent node isTip
#> 51 139 51 TRUE
#> 52 139 52 TRUE
#> 63 147 63 TRUE
#> 67 150 67 TRUE
#> 136 135 136 FALSE
#> 71 150 71 TRUE
#> 125 100 125 FALSE
#> 20 115 20 TRUE
#> 116 115 116 FALSE
#> 140 100 140 FALSE
#> clusters
#> 51 77
#> 52 66
#> 63 41
#> 67 32
#> 136 58, 67, 57
#> 71 54
#> 125 85, 84, 95, 97, 87, 74, 76, 96, 86, 69, 78, 56, 77, 66, 75, 68, 58, 67, 57
#> 20 89
#> 116 59, 88, 70, 60
#> 140 35, 43, 46, 73, 64, 100, 45, 21, 42, 31, 41, 33, 81, 34, 32, 52, 55, 44, 54, 53, 63, 24, 82, 23, 61, 51, 71, 62
#> stat_total stat_parent pval_total pval_parent
#> 51 1.0823881 4.439944 0.290116635 0.0003193122
#> 52 -1.9247385 -4.439944 0.069760334 0.0003193122
#> 63 1.8107720 3.337019 0.080505219 0.0022734266
#> 67 1.7508315 3.229075 0.092603817 0.0033277959
#> 136 1.4290582 -3.198607 0.163963670 0.0041906737
#> 71 -2.3514750 -2.829766 0.033743471 0.0082310590
#> 125 2.7073249 2.707325 0.011833319 0.0118333189
#> 20 -0.7764417 -2.897438 0.446585533 0.0124821272
#> 116 3.0537776 2.897438 0.005462897 0.0124821272
#> 140 -2.5102007 -2.510201 0.017809639 0.0178096392The results of the previous steps are visualised by a coloured tree with a corresponding heatmap. The heatmap displays the median scaled marker expressions of each cluster to help understand what cell type each cluster may represent, and the tree not only reveals how clusters have been hierarchically aggregated, but is coloured on each node by the test statistic obtained when testing using the proportions relative to all of that node, with the branch connecting the child to the parent coloured by the test statistic obtained when testing using the proportions relative to parent of the child node.
plotInteractiveHeatmap(tested_tree,
clust_med_df = clust_tree$median_freq,
clusters=clusters)
#> Warning: Removed 1 rows containing missing values (geom_interactive_point).Below we change the hierarchical aggregation technique to average-linkage hierarchical clustering.
The available options include any of the available ethods in the hclust() function, ie, one of "ward.D", "ward.D2", "single", "complete", "average", "mcquitty", "median" or "centroid"
clust_tree <- getClusterTree(exprs,
clusters,
hierarchy_method="average")
tested_tree <- testTree(clust_tree$clust_tree,
clusters=clusters,
samples=samples,
classes=classes,
pos_class_name=NULL)
plotInteractiveHeatmap(tested_tree,
clust_med_df = clust_tree$median_freq,
clusters=clusters)
#> Warning: Removed 1 rows containing missing values (geom_interactive_point).Below we change the significance test to use count models instead of a t-test/Wilcoxon test. The available methods are "GLMM" or "edgeR"
clust_tree <- getClusterTree(exprs,
clusters,
hierarchy_method="hopach")
tested_tree_edgeR <- testTree(clust_tree$clust_tree,
clusters=clusters,
samples=samples,
classes=classes,
sig_test="edgeR",
pos_class_name="COV")
#> Using classic mode.
plotInteractiveHeatmap(tested_tree_edgeR,
clust_med_df = clust_tree$median_freq,
clusters=clusters)You can also extract edgeR FDR values as so:
head(tested_tree_edgeR$data[,c("parent", "node", "isTip", "clusters", "FDR_parent", "FDR_total")])
#> # A tibble: 6 × 6
#> parent node isTip clusters FDR_parent FDR_total
#> <int> <dbl> <lgl> <list> <dbl> <dbl>
#> 1 102 1 TRUE <chr [1]> 0.890 0.158
#> 2 102 2 TRUE <chr [1]> 0.808 0.157
#> 3 104 3 TRUE <chr [1]> 0.794 0.920
#> 4 104 4 TRUE <chr [1]> 0.278 0.107
#> 5 105 5 TRUE <chr [1]> 0.794 0.331
#> 6 105 6 TRUE <chr [1]> 0.943 0.636Extract cell type proportions (both %total and %parent), as well as the geometric mean of each marker per cell type. These can then be used for further investigation via visualisation or machine learning.
The function below returns a dataframe containing the absolute proportions and proportions to parent for each cell type for each sample.
perc_parent_n_celltype denotes the proportion of the cell type relative to its parent (%parent) n generations away, ie. perc_parent_1_2 is the proportion of cluster 2 relative to its direct parent in the hierarchical tree.
perc_total_celltype denotes the proportion of the cell type relative all cells (%total)
prop_df <- getCellProp(phylo=clust_tree$clust_tree,
clusters=clusters,
samples=samples,
classes=classes)
head(prop_df[,1:8])
#> sample_id class perc_total_92 perc_total_91 perc_total_7 perc_total_50
#> 1 COV1 COV 0.029411765 0.005882353 0.01176471 0.000000000
#> 2 COV10 COV 0.066666667 0.044444444 0.00000000 0.005555556
#> 3 COV12 COV 0.051428571 0.000000000 0.03428571 0.000000000
#> 4 COV13 COV 0.006410256 0.000000000 0.01923077 0.012820513
#> 5 COV14 COV 0.019867550 0.019867550 0.00000000 0.006622517
#> 6 COV15 COV 0.000000000 0.000000000 0.00000000 0.000000000
#> perc_total_5 perc_total_39
#> 1 0.005882353 0.000000000
#> 2 0.016666667 0.005555556
#> 3 0.005714286 0.000000000
#> 4 0.038461538 0.000000000
#> 5 0.006622517 0.000000000
#> 6 0.006993007 0.000000000The function below returns a dataframe containing the geometric mean for each marker for each cell type for each sample.
m_gmean_celltype denotes the geometric mean of marker m in the cell type per sample
means_df <- getCellGMeans(clust_tree$clust_tree,
exprs=exprs,
clusters=clusters,
samples=samples,
classes=classes)
head(means_df[,1:8])
#> sample_id class CD45RA B525-FITC-H_gmean_92 CD45RA B525-FITC-H_gmean_91
#> 1 COV1 COV 8.026009 10.271191
#> 2 COV10 COV 8.804902 8.432155
#> 3 COV12 COV 8.268734 NaN
#> 4 COV13 COV 10.221825 NaN
#> 5 COV14 COV 9.768303 10.470112
#> 6 COV15 COV NaN NaN
#> CD45RA B525-FITC-H_gmean_7 CD45RA B525-FITC-H_gmean_50
#> 1 7.551528 NaN
#> 2 NaN 6.360800
#> 3 6.611482 NaN
#> 4 8.369937 10.629523
#> 5 NaN 8.597188
#> 6 NaN NaN
#> CD45RA B525-FITC-H_gmean_5 CD45RA B525-FITC-H_gmean_39
#> 1 7.126569 NaN
#> 2 8.496825 7.683541
#> 3 9.265149 NaN
#> 4 8.178744 NaN
#> 5 8.543981 NaN
#> 6 11.036683 NaNsessionInfo()
#> R version 4.1.2 (2021-11-01)
#> Platform: x86_64-pc-linux-gnu (64-bit)
#> Running under: Ubuntu 20.04.3 LTS
#>
#> Matrix products: default
#> BLAS: /home/biocbuild/bbs-3.14-bioc/R/lib/libRblas.so
#> LAPACK: /home/biocbuild/bbs-3.14-bioc/R/lib/libRlapack.so
#>
#> 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
#>
#> attached base packages:
#> [1] stats4 stats graphics grDevices utils datasets methods
#> [8] base
#>
#> other attached packages:
#> [1] ggtree_3.2.1 SingleCellExperiment_1.16.0
#> [3] SummarizedExperiment_1.24.0 Biobase_2.54.0
#> [5] GenomicRanges_1.46.1 GenomeInfoDb_1.30.1
#> [7] IRanges_2.28.0 S4Vectors_0.32.3
#> [9] BiocGenerics_0.40.0 MatrixGenerics_1.6.0
#> [11] matrixStats_0.61.0 treekoR_1.2.1
#> [13] BiocStyle_2.22.0
#>
#> loaded via a namespace (and not attached):
#> [1] utf8_1.2.2 tidyselect_1.1.1
#> [3] lme4_1.1-28 htmlwidgets_1.5.4
#> [5] BiocParallel_1.28.3 grid_4.1.2
#> [7] Rtsne_0.15 aws.signature_0.6.0
#> [9] flowCore_2.6.0 ScaledMatrix_1.2.0
#> [11] munsell_0.5.0 codetools_0.2-18
#> [13] colorspace_2.0-2 knitr_1.37
#> [15] uuid_1.0-3 ggsignif_0.6.3
#> [17] labeling_0.4.2 GenomeInfoDbData_1.2.7
#> [19] polyclip_1.10-0 farver_2.1.0
#> [21] pheatmap_1.0.12 flowWorkspace_4.6.0
#> [23] vctrs_0.3.8 treeio_1.18.1
#> [25] generics_0.1.2 TH.data_1.1-0
#> [27] xfun_0.29 R6_2.5.1
#> [29] doParallel_1.0.17 ggbeeswarm_0.6.0
#> [31] clue_0.3-60 rsvd_1.0.5
#> [33] locfit_1.5-9.4 ggiraph_0.8.1
#> [35] bitops_1.0-7 gridGraphics_0.5-1
#> [37] DelayedArray_0.20.0 assertthat_0.2.1
#> [39] scales_1.1.1 multcomp_1.4-18
#> [41] beeswarm_0.4.0 gtable_0.3.0
#> [43] beachmat_2.10.0 RProtoBufLib_2.6.0
#> [45] sandwich_3.0-1 rlang_1.0.1
#> [47] systemfonts_1.0.4 GlobalOptions_0.1.2
#> [49] splines_4.1.2 rstatix_0.7.0
#> [51] lazyeval_0.2.2 hexbin_1.28.2
#> [53] broom_0.7.12 BiocManager_1.30.16
#> [55] yaml_2.2.2 reshape2_1.4.4
#> [57] abind_1.4-5 backports_1.4.1
#> [59] RBGL_1.70.0 tools_4.1.2
#> [61] bookdown_0.24 ggplotify_0.1.0
#> [63] ggplot2_3.3.5 ellipsis_0.3.2
#> [65] jquerylib_0.1.4 RColorBrewer_1.1-2
#> [67] ggridges_0.5.3 Rcpp_1.0.8
#> [69] plyr_1.8.6 sparseMatrixStats_1.6.0
#> [71] base64enc_0.1-3 zlibbioc_1.40.0
#> [73] purrr_0.3.4 RCurl_1.98-1.6
#> [75] FlowSOM_2.2.0 ggpubr_0.4.0
#> [77] viridis_0.6.2 GetoptLong_1.0.5
#> [79] cowplot_1.1.1 zoo_1.8-9
#> [81] CATALYST_1.18.1 ggrepel_0.9.1
#> [83] cluster_2.1.2 colorRamps_2.3
#> [85] magrittr_2.0.2 ncdfFlow_2.40.0
#> [87] data.table_1.14.2 scattermore_0.7
#> [89] circlize_0.4.14 mvtnorm_1.1-3
#> [91] ggnewscale_0.4.5 patchwork_1.1.1
#> [93] evaluate_0.14 XML_3.99-0.8
#> [95] jpeg_0.1-9 gridExtra_2.3
#> [97] shape_1.4.6 ggcyto_1.22.0
#> [99] scater_1.22.0 compiler_4.1.2
#> [101] tibble_3.1.6 crayon_1.4.2
#> [103] ggpointdensity_0.1.0 minqa_1.2.4
#> [105] htmltools_0.5.2 ggfun_0.0.5
#> [107] tidyr_1.2.0 aplot_0.1.2
#> [109] RcppParallel_5.1.5 aws.s3_0.3.21
#> [111] DBI_1.1.2 tweenr_1.0.2
#> [113] ComplexHeatmap_2.10.0 MASS_7.3-55
#> [115] boot_1.3-28 Matrix_1.4-0
#> [117] car_3.0-12 diffcyt_1.14.0
#> [119] cli_3.1.1 parallel_4.1.2
#> [121] igraph_1.2.11 pkgconfig_2.0.3
#> [123] scuttle_1.4.0 hopach_2.54.0
#> [125] xml2_1.3.3 foreach_1.5.2
#> [127] vipor_0.4.5 bslib_0.3.1
#> [129] XVector_0.34.0 drc_3.0-1
#> [131] yulab.utils_0.0.4 stringr_1.4.0
#> [133] digest_0.6.29 ConsensusClusterPlus_1.58.0
#> [135] graph_1.72.0 rmarkdown_2.11
#> [137] tidytree_0.3.7 edgeR_3.36.0
#> [139] DelayedMatrixStats_1.16.0 curl_4.3.2
#> [141] gtools_3.9.2 rjson_0.2.21
#> [143] nloptr_2.0.0 lifecycle_1.0.1
#> [145] nlme_3.1-155 jsonlite_1.7.3
#> [147] BiocNeighbors_1.12.0 carData_3.0-5
#> [149] viridisLite_0.4.0 limma_3.50.0
#> [151] fansi_1.0.2 pillar_1.7.0
#> [153] lattice_0.20-45 plotrix_3.8-2
#> [155] fastmap_1.1.0 httr_1.4.2
#> [157] survival_3.2-13 glue_1.6.1
#> [159] png_0.1-7 iterators_1.0.14
#> [161] Rgraphviz_2.38.0 nnls_1.4
#> [163] ggforce_0.3.3 stringi_1.7.6
#> [165] sass_0.4.0 BiocSingular_1.10.0
#> [167] CytoML_2.6.0 latticeExtra_0.6-29
#> [169] dplyr_1.0.8 cytolib_2.6.2
#> [171] irlba_2.3.5 ape_5.6-1