# treekoR

Adam Chan1* and Ellis Patrick1,2

1School of Mathematics and Statistics, The University of Sydney, Sydney, New South Wales, Australia
2Westmead Institute for Medical Research, University of Sydney, Australia

*adam.s.chan@sydney.edu.au

#### 15 February 2022

# Contents

    # Load required packages
    suppressPackageStartupMessages({
      library(treekoR)
      library(SingleCellExperiment)
      library(ggtree)
    })

# 1 Installation

    # Install the development version from GitHub:
    # install.packages("devtools")
    devtools::install_github("adam2o1o/treekoR")
    library(treekoR)

# 2 Overview

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.

# 3 Usage

## 3.1 Example Data

- DeBiasi_COVID_CD8_samp: A `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_samp

## 3.2 Set up

treekoR 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 cells
- `clusters`: 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 to
- `samples`: a vector of length \(n\) identifying the patient each cell belongs to

In 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_id

## 3.3 Construction of Hierarchy of Clusters

The 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")

## 3.4 Significance testing of Cell Subpopulations

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)

### 3.4.1 Extracting significance results

- `node`: unique identifier for each node in the hierarchical tree
- `parent`: the node of the parent
- `isTip`: whether the node is a leaf node in the tree
- `clusters`: the clusters belonging to the corresponding node
- `stat_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.0178096392

## 3.5 Interactive Visualisation

The 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).

### 3.5.1 Using different hierarchical aggregations

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).

### 3.5.2 Using counts models for significance testing

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 x 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.636

## 3.6 Feature extraction

Extract 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.

### 3.6.1 Extracting Proportions

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.000000000

### 3.6.2 Extracting Expression Geometric Means

The 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                         NaN

    sessionInfo()
    #> R version 4.1.2 (2021-11-01)
    #> Platform: x86_64-w64-mingw32/x64 (64-bit)
    #> Running under: Windows Server 2012 R2 x64 (build 9600)
    #>
    #> Matrix products: default
    #>
    #> locale:
    #> [1] LC_COLLATE=C
    #> [2] LC_CTYPE=English_United States.1252
    #> [3] LC_MONETARY=English_United States.1252
    #> [4] LC_NUMERIC=C
    #> [5] LC_TIME=English_United States.1252
    #>
    #> 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.8
    #>  [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.5.0
    #> [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.2.0                   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