Feature extraction
Timothy Keyes
2024-05-07
Source:vignettes/feature-extraction.Rmd
feature-extraction.Rmd
library(tidytof)
library(dplyr)
#>
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:stats':
#>
#> filter, lag
#> The following objects are masked from 'package:base':
#>
#> intersect, setdiff, setequal, union
library(stringr)In addition to its functions for analyzing and visualizing CyTOF data
at the single-cell and cluster levels, tidytof’s
tof_extract_features() verb allows users to aggregate
single-cell and cluster-level information in order to summarize whole
samples (or whole patients) from which cells were collected. These
features can be useful for visualizing the differences between patients
and samples in different experimental conditions or for building machine
learning models.
To understand how the tof_extract_features() verb works,
it’s easiest to look at each of its subroutines (the members of the
tof_extract_* function family) independently.
Accessing the data for this vignette
To demonstrate how to use these verbs, we’ll first download a dataset originally collected for the development of the CITRUS algorithm. These data are available in the HDCytoData package, which is available on Bioconductor and can be downloaded with the following command:
if (!requireNamespace("BiocManager", quietly = TRUE)) {
install.packages("BiocManager")
}
BiocManager::install("HDCytoData")To load the CITRUS data into our current R session, we can call a
function from the HDCytoData, which will provide it to us
in a format from the {flowCore} package (called a
“flowSet”). To convert this into a tidy tibble, we can use
tidytof built-in method for converting flowCore objects
into tof_tbl’s .
citrus_raw <- HDCytoData::Bodenmiller_BCR_XL_flowSet()
citrus_data <-
citrus_raw |>
as_tof_tbl(sep = "_")Thus, we can see that citrus_data is a
tof_tbl with 172791 cells (one in each row) and 39 pieces
of information about each cell (one in each column).
We can also extract some metadata from the raw data and join it with
our single-cell data using some functions from the
tidyverse:
citrus_metadata <-
tibble(
file_name = as.character(flowCore::pData(citrus_raw)[[1]]),
sample_id = 1:length(file_name),
patient = stringr::str_extract(file_name, "patient[:digit:]"),
stimulation = stringr::str_extract(file_name, "(BCR-XL)|Reference")
) |>
mutate(
stimulation = if_else(stimulation == "Reference", "Basal", stimulation)
)
citrus_metadata |>
head()
#> # A tibble: 6 × 4
#> file_name sample_id patient stimulation
#> <chr> <int> <chr> <chr>
#> 1 PBMC8_30min_patient1_BCR-XL.fcs 1 patient1 BCR-XL
#> 2 PBMC8_30min_patient1_Reference.fcs 2 patient1 Basal
#> 3 PBMC8_30min_patient2_BCR-XL.fcs 3 patient2 BCR-XL
#> 4 PBMC8_30min_patient2_Reference.fcs 4 patient2 Basal
#> 5 PBMC8_30min_patient3_BCR-XL.fcs 5 patient3 BCR-XL
#> 6 PBMC8_30min_patient3_Reference.fcs 6 patient3 BasalThus, we now have sample-level information about which patient each sample was collected from and which stimulation condition (“Basal” or “BCR-XL”) each sample was exposed to before data acquisition.
Finally, we can join this metadata with our single-cell
tof_tbl to obtain the cleaned dataset.
citrus_data <-
citrus_data |>
left_join(citrus_metadata, by = "sample_id")After these data cleaning steps, we now have
citrus_data, a tof_tbl containing cells
collected from 8 patients. Specifically, 2 samples were taken from each
patient: one in which the cells’ B-cell receptors were stimulated
(BCR-XL) and one in which they were not (Basal). In
citrus_data, each cell’s patient of origin is stored in the
patient column, and each cell’s stimulation condition is
stored in the stimulation column. In addition, the
population_id column stores information about cluster
labels that were applied to each cell using a combination of FlowSOM
clustering and manual merging (for details, run
?HDCytoData::Bodenmiller_BCR_XL in the R console).
Calculating cluster proportions using
tof_extract_proportion()
First, we have tof_extract_proportion(), which extracts
the proportion of cells in each cluster within each sample (with samples
defined using the group_cols argument):
# preprocess the numeric columns in the citrus dataset
citrus_data <-
citrus_data |>
mutate(cluster = str_c("cluster", population_id)) |>
tof_preprocess()
citrus_data |>
tof_extract_proportion(
cluster_col = cluster,
group_cols = c(patient, stimulation)
) |>
head()
#> # A tibble: 6 × 10
#> patient stimulation `prop@cluster1` `prop@cluster2` `prop@cluster3`
#> <chr> <chr> <dbl> <dbl> <dbl>
#> 1 patient1 Basal 0.0190 0.0482 0.447
#> 2 patient1 BCR-XL 0.0109 0.0395 0.268
#> 3 patient2 Basal 0.0130 0.0280 0.491
#> 4 patient2 BCR-XL 0.0101 0.0143 0.358
#> 5 patient3 Basal 0.0326 0.0830 0.397
#> 6 patient3 BCR-XL 0.0200 0.0412 0.323
#> # ℹ 5 more variables: `prop@cluster4` <dbl>, `prop@cluster5` <dbl>,
#> # `prop@cluster6` <dbl>, `prop@cluster7` <dbl>, `prop@cluster8` <dbl>Like all members of the tof_extract_* function family,
tof_extract_proportion() returns one row for each sample
(defined as a unique combination of values of the columns specified in
group_cols) and one column for each extracted feature
(above, one column for the proportion of each of the 8 clusters in
citrus_data). These values can also be returned in “long”
format by changing the format argument:
citrus_data |>
tof_extract_proportion(
cluster_col = cluster,
group_cols = c(patient, stimulation),
format = "long"
) |>
head()
#> # A tibble: 6 × 4
#> patient stimulation cluster prop
#> <chr> <chr> <chr> <dbl>
#> 1 patient1 Basal cluster1 0.0190
#> 2 patient1 Basal cluster2 0.0482
#> 3 patient1 Basal cluster3 0.447
#> 4 patient1 Basal cluster4 0.237
#> 5 patient1 Basal cluster5 0.00219
#> 6 patient1 Basal cluster6 0.0759Calculating cluster marker expression measures using
tof_extract_central_tendency()
Another member of the tof_extract_*() function family,
tof_extract_central_tendency(), computes the central
tendency (e.g. mean or median) of user-specified markers in each
cluster.
citrus_data |>
tof_extract_central_tendency(
cluster_col = cluster,
group_cols = c(patient, stimulation),
marker_cols = any_of(c("CD45_In115", "CD4_Nd145", "CD20_Sm147")),
central_tendency_function = mean
) |>
head()
#> # A tibble: 6 × 26
#> patient stimulation `CD45_In115@cluster1_ct` `CD4_Nd145@cluster1_ct`
#> <chr> <chr> <dbl> <dbl>
#> 1 patient1 BCR-XL 4.80 0.0967
#> 2 patient1 Basal 4.68 0.765
#> 3 patient2 BCR-XL 5.00 -0.0579
#> 4 patient2 Basal 4.88 0.808
#> 5 patient3 BCR-XL 5.04 -0.0432
#> 6 patient3 Basal 4.98 0.745
#> # ℹ 22 more variables: `CD20_Sm147@cluster1_ct` <dbl>,
#> # `CD45_In115@cluster2_ct` <dbl>, `CD4_Nd145@cluster2_ct` <dbl>,
#> # `CD20_Sm147@cluster2_ct` <dbl>, `CD45_In115@cluster3_ct` <dbl>,
#> # `CD4_Nd145@cluster3_ct` <dbl>, `CD20_Sm147@cluster3_ct` <dbl>,
#> # `CD45_In115@cluster4_ct` <dbl>, `CD4_Nd145@cluster4_ct` <dbl>,
#> # `CD20_Sm147@cluster4_ct` <dbl>, `CD45_In115@cluster5_ct` <dbl>,
#> # `CD4_Nd145@cluster5_ct` <dbl>, `CD20_Sm147@cluster5_ct` <dbl>, …The argument central_tendency_function can be used to
compute any summary statistic. For example, the following choice for
central_tendency_function will compute the 75th percentile
for each marker-cluster pair in citrus_data:
citrus_data |>
tof_extract_central_tendency(
cluster_col = cluster,
group_cols = c(patient, stimulation),
marker_cols = any_of(c("CD45_In115", "CD4_Nd145", "CD20_Sm147")),
central_tendency_function = function(x) quantile(x = x, probs = 0.75)
) |>
head()
#> # A tibble: 6 × 26
#> patient stimulation `CD45_In115@cluster1_ct` `CD4_Nd145@cluster1_ct`
#> <chr> <chr> <dbl> <dbl>
#> 1 patient1 BCR-XL 5.30 -0.0186
#> 2 patient1 Basal 5.18 1.32
#> 3 patient2 BCR-XL 5.41 -0.0201
#> 4 patient2 Basal 5.28 1.39
#> 5 patient3 BCR-XL 5.42 -0.0362
#> 6 patient3 Basal 5.41 1.27
#> # ℹ 22 more variables: `CD20_Sm147@cluster1_ct` <dbl>,
#> # `CD45_In115@cluster2_ct` <dbl>, `CD4_Nd145@cluster2_ct` <dbl>,
#> # `CD20_Sm147@cluster2_ct` <dbl>, `CD45_In115@cluster3_ct` <dbl>,
#> # `CD4_Nd145@cluster3_ct` <dbl>, `CD20_Sm147@cluster3_ct` <dbl>,
#> # `CD45_In115@cluster4_ct` <dbl>, `CD4_Nd145@cluster4_ct` <dbl>,
#> # `CD20_Sm147@cluster4_ct` <dbl>, `CD45_In115@cluster5_ct` <dbl>,
#> # `CD4_Nd145@cluster5_ct` <dbl>, `CD20_Sm147@cluster5_ct` <dbl>, …Calculating the proportion of cells with marker expression above a
threshold using tof_extract_proportion()
tof_extract_threshold() is similar to
tof_extract_central_tendency(), but calculates the
proportion of cells above a user-specified expression value for each
marker instead of a measure of central tendency:
citrus_data |>
tof_extract_threshold(
cluster_col = cluster,
group_cols = c(patient, stimulation),
marker_cols = any_of(c("CD45_In115", "CD4_Nd145", "CD20_Sm147")),
threshold = 5
) |>
head()
#> # A tibble: 6 × 26
#> patient stimulation `CD45_In115@cluster1_threshold` CD4_Nd145@cluster1_thre…¹
#> <chr> <chr> <dbl> <dbl>
#> 1 patient1 BCR-XL 0.516 0
#> 2 patient1 Basal 0.365 0
#> 3 patient2 BCR-XL 0.554 0
#> 4 patient2 Basal 0.452 0
#> 5 patient3 BCR-XL 0.547 0
#> 6 patient3 Basal 0.549 0
#> # ℹ abbreviated name: ¹`CD4_Nd145@cluster1_threshold`
#> # ℹ 22 more variables: `CD20_Sm147@cluster1_threshold` <dbl>,
#> # `CD45_In115@cluster2_threshold` <dbl>,
#> # `CD4_Nd145@cluster2_threshold` <dbl>,
#> # `CD20_Sm147@cluster2_threshold` <dbl>,
#> # `CD45_In115@cluster3_threshold` <dbl>,
#> # `CD4_Nd145@cluster3_threshold` <dbl>, …Calculating differences in marker distributions using
tof_extract_emd() and tof_extract_jsd()
The two final members of the tof_extract_* function
family – tof_extract_emd and tof_extract_jsd –
are designed specifically for comparing distributions of marker
expression between stimulation conditions. As such, they must be given a
stimulation column (using the emd_col or
jsd_col argument) that identifies the stimulation condition
each cell is in, and a reference_level that specifies the
reference (i.e. unstimulated) condition within the emd_col
or jsd_col.
With these additional arguments, tof_extract_emd
computes the Earth-mover’s distance between each marker’s distribution
in the stimulation conditions (within each cluster) and the basal
condition; similarly, tof_extract_jsd computes the
Jensen-Shannon divergence index between the same distributions. Both of
these values are ways to compare how different 2 distributions are to
one another and are more computationally expensive (but also
higher-resolution) than simply comparing measures of central
tendency.
# Earth-mover's distance
citrus_data |>
tof_extract_emd(
cluster_col = cluster,
group_cols = patient,
marker_cols = any_of(c("CD45_In115", "CD4_Nd145", "CD20_Sm147")),
emd_col = stimulation,
reference_level = "Basal"
) |>
head()
#> # A tibble: 6 × 25
#> patient BCR-XL_CD45_In115@clu…¹ BCR-XL_CD4_Nd145@clu…² BCR-XL_CD20_Sm147@cl…³
#> <chr> <dbl> <dbl> <dbl>
#> 1 patient1 0.864 2.47 13.0
#> 2 patient2 1.11 7.05 10.8
#> 3 patient3 0.670 6.23 10.5
#> 4 patient4 2.64 5.86 9.90
#> 5 patient5 0.594 7.56 8.13
#> 6 patient6 0.661 4.77 7.97
#> # ℹ abbreviated names: ¹`BCR-XL_CD45_In115@cluster3_emd`,
#> # ²`BCR-XL_CD4_Nd145@cluster3_emd`, ³`BCR-XL_CD20_Sm147@cluster3_emd`
#> # ℹ 21 more variables: `BCR-XL_CD45_In115@cluster7_emd` <dbl>,
#> # `BCR-XL_CD4_Nd145@cluster7_emd` <dbl>,
#> # `BCR-XL_CD20_Sm147@cluster7_emd` <dbl>,
#> # `BCR-XL_CD45_In115@cluster4_emd` <dbl>,
#> # `BCR-XL_CD4_Nd145@cluster4_emd` <dbl>, …
# Jensen-Shannon Divergence
citrus_data |>
tof_extract_jsd(
cluster_col = cluster,
group_cols = patient,
marker_cols = any_of(c("CD45_In115", "CD4_Nd145", "CD20_Sm147")),
jsd_col = stimulation,
reference_level = "Basal"
) |>
head()
#> # A tibble: 6 × 25
#> patient BCR-XL_CD45_In115@clu…¹ BCR-XL_CD4_Nd145@clu…² BCR-XL_CD20_Sm147@cl…³
#> <chr> <dbl> <dbl> <dbl>
#> 1 patient1 0.0367 0.0513 0.347
#> 2 patient2 0.00831 0.168 0.401
#> 3 patient3 0.0104 0.115 0.357
#> 4 patient4 0.0301 0.135 0.205
#> 5 patient5 0.00911 0.0789 0.274
#> 6 patient6 0.00972 0.0346 0.214
#> # ℹ abbreviated names: ¹`BCR-XL_CD45_In115@cluster3_jsd`,
#> # ²`BCR-XL_CD4_Nd145@cluster3_jsd`, ³`BCR-XL_CD20_Sm147@cluster3_jsd`
#> # ℹ 21 more variables: `BCR-XL_CD45_In115@cluster7_jsd` <dbl>,
#> # `BCR-XL_CD4_Nd145@cluster7_jsd` <dbl>,
#> # `BCR-XL_CD20_Sm147@cluster7_jsd` <dbl>,
#> # `BCR-XL_CD45_In115@cluster4_jsd` <dbl>,
#> # `BCR-XL_CD4_Nd145@cluster4_jsd` <dbl>, …Putting it all together with
tof_extract_features()
Finally, the tof_extract_features() verb provides a
wrapper to each of the members of its function family, allowing users to
extract multiple features types at once. For example, the following code
extracts the proportion of each cluster, median of several markers in
each cluster, and EMD between the basal condition and stimulated
condition in each cluster for all patients in
citrus_data.
signaling_markers <-
c(
"pNFkB_Nd142", "pStat5_Nd150", "pAkt_Sm152", "pStat1_Eu153", "pStat3_Gd158",
"pSlp76_Dy164", "pBtk_Er166", "pErk_Er168", "pS6_Yb172", "pZap70_Gd156"
)
citrus_data |>
tof_extract_features(
cluster_col = cluster,
group_cols = patient,
stimulation_col = stimulation,
lineage_cols = any_of(c("CD45_In115", "CD20_Sm147", "CD33_Nd148")),
signaling_cols = any_of(signaling_markers),
signaling_method = "emd",
basal_level = "Basal"
) |>
head()
#> # A tibble: 6 × 193
#> patient `prop@cluster1` `prop@cluster2` `prop@cluster3` `prop@cluster4`
#> <chr> <dbl> <dbl> <dbl> <dbl>
#> 1 patient1 0.0149 0.0438 0.356 0.351
#> 2 patient2 0.0115 0.0212 0.425 0.323
#> 3 patient3 0.0255 0.0594 0.355 0.217
#> 4 patient4 0.0127 0.0418 0.320 0.223
#> 5 patient5 0.0207 0.0423 0.377 0.269
#> 6 patient6 0.0183 0.0493 0.459 0.250
#> # ℹ 188 more variables: `prop@cluster5` <dbl>, `prop@cluster6` <dbl>,
#> # `prop@cluster7` <dbl>, `prop@cluster8` <dbl>,
#> # `CD45_In115@cluster1_ct` <dbl>, `CD20_Sm147@cluster1_ct` <dbl>,
#> # `CD33_Nd148@cluster1_ct` <dbl>, `CD45_In115@cluster2_ct` <dbl>,
#> # `CD20_Sm147@cluster2_ct` <dbl>, `CD33_Nd148@cluster2_ct` <dbl>,
#> # `CD45_In115@cluster3_ct` <dbl>, `CD20_Sm147@cluster3_ct` <dbl>,
#> # `CD33_Nd148@cluster3_ct` <dbl>, `CD45_In115@cluster4_ct` <dbl>, …Session info
sessionInfo()
#> R version 4.4.0 (2024-04-24)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 22.04.4 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.20.so; LAPACK version 3.10.0
#>
#> locale:
#> [1] LC_CTYPE=C.UTF-8 LC_NUMERIC=C LC_TIME=C.UTF-8
#> [4] LC_COLLATE=C.UTF-8 LC_MONETARY=C.UTF-8 LC_MESSAGES=C.UTF-8
#> [7] LC_PAPER=C.UTF-8 LC_NAME=C LC_ADDRESS=C
#> [10] LC_TELEPHONE=C LC_MEASUREMENT=C.UTF-8 LC_IDENTIFICATION=C
#>
#> time zone: UTC
#> tzcode source: system (glibc)
#>
#> attached base packages:
#> [1] stats4 stats graphics grDevices utils datasets methods
#> [8] base
#>
#> other attached packages:
#> [1] HDCytoData_1.24.0 flowCore_2.16.0
#> [3] SummarizedExperiment_1.34.0 Biobase_2.64.0
#> [5] GenomicRanges_1.56.0 GenomeInfoDb_1.40.0
#> [7] IRanges_2.38.0 S4Vectors_0.42.0
#> [9] MatrixGenerics_1.16.0 matrixStats_1.3.0
#> [11] ExperimentHub_2.12.0 AnnotationHub_3.12.0
#> [13] BiocFileCache_2.12.0 dbplyr_2.5.0
#> [15] BiocGenerics_0.50.0 stringr_1.5.1
#> [17] dplyr_1.1.4 tidytof_0.99.7
#>
#> loaded via a namespace (and not attached):
#> [1] jsonlite_1.8.8 shape_1.4.6.1 magrittr_2.0.3
#> [4] farver_2.1.1 rmarkdown_2.26 fs_1.6.4
#> [7] zlibbioc_1.50.0 ragg_1.3.0 vctrs_0.6.5
#> [10] memoise_2.0.1 htmltools_0.5.8.1 S4Arrays_1.4.0
#> [13] curl_5.2.1 SparseArray_1.4.0 sass_0.4.9
#> [16] parallelly_1.37.1 bslib_0.7.0 htmlwidgets_1.6.4
#> [19] desc_1.4.3 lubridate_1.9.3 cachem_1.0.8
#> [22] igraph_2.0.3 mime_0.12 lifecycle_1.0.4
#> [25] iterators_1.0.14 pkgconfig_2.0.3 Matrix_1.7-0
#> [28] R6_2.5.1 fastmap_1.1.1 GenomeInfoDbData_1.2.12
#> [31] future_1.33.2 digest_0.6.35 colorspace_2.1-0
#> [34] AnnotationDbi_1.66.0 textshaping_0.3.7 RSQLite_2.3.6
#> [37] philentropy_0.8.0 filelock_1.0.3 cytolib_2.16.0
#> [40] fansi_1.0.6 yardstick_1.3.1 timechange_0.3.0
#> [43] httr_1.4.7 polyclip_1.10-6 abind_1.4-5
#> [46] compiler_4.4.0 bit64_4.0.5 withr_3.0.0
#> [49] doParallel_1.0.17 viridis_0.6.5 DBI_1.2.2
#> [52] ggforce_0.4.2 MASS_7.3-60.2 lava_1.8.0
#> [55] rappdirs_0.3.3 DelayedArray_0.30.0 tools_4.4.0
#> [58] future.apply_1.11.2 nnet_7.3-19 glue_1.7.0
#> [61] grid_4.4.0 generics_0.1.3 recipes_1.0.10
#> [64] gtable_0.3.5 tzdb_0.4.0 class_7.3-22
#> [67] tidyr_1.3.1 data.table_1.15.4 hms_1.1.3
#> [70] tidygraph_1.3.1 utf8_1.2.4 XVector_0.44.0
#> [73] ggrepel_0.9.5 BiocVersion_3.19.1 foreach_1.5.2
#> [76] pillar_1.9.0 RcppHNSW_0.6.0 splines_4.4.0
#> [79] tweenr_2.0.3 lattice_0.22-6 survival_3.5-8
#> [82] bit_4.0.5 emdist_0.3-3 RProtoBufLib_2.16.0
#> [85] tidyselect_1.2.1 Biostrings_2.72.0 knitr_1.46
#> [88] gridExtra_2.3 xfun_0.43 graphlayouts_1.1.1
#> [91] hardhat_1.3.1 timeDate_4032.109 stringi_1.8.3
#> [94] UCSC.utils_1.0.0 yaml_2.3.8 evaluate_0.23
#> [97] codetools_0.2-20 ggraph_2.2.1 tibble_3.2.1
#> [100] BiocManager_1.30.23 cli_3.6.2 rpart_4.1.23
#> [103] systemfonts_1.0.6 munsell_0.5.1 jquerylib_0.1.4
#> [106] Rcpp_1.0.12 globals_0.16.3 png_0.1-8
#> [109] parallel_4.4.0 pkgdown_2.0.9 gower_1.0.1
#> [112] ggplot2_3.5.1 readr_2.1.5 blob_1.2.4
#> [115] listenv_0.9.1 glmnet_4.1-8 viridisLite_0.4.2
#> [118] ipred_0.9-14 scales_1.3.0 prodlim_2023.08.28
#> [121] purrr_1.0.2 crayon_1.5.2 rlang_1.1.3
#> [124] KEGGREST_1.44.0