Censored covariate in cytometry
29 March 2022
Abstract
censcyt provides functionalities for differential abundance analysis in cytometry when a covariate is subject to right censoring (e.g. survival time).
Package
censcyt 1.2.1
Contents
1 Introduction
The censcyt package provides functionalities for differential abundance analysis in cytometry when a covariate is subject to right censoring. It extends the differential
discovery capabilities of the diffcyt package (Weber et al. 2019)
by including association testing between the cell population abundance and a censored variable such as survival time. The general workflow is the same as described in the diffcyt package vignette
and it is advisable to be familiar with it before continuing.
As a short summary, the diffcyt workflow consists of preprocessing of raw cytometry data, cell type assignment through automated clustering and differential abundance analysis (or differential state analysis) (Weber et al. 2019).
If a time-to-event variable (e.g. survival time) connected to a cytometry sample is available, one might be interested in the association between the cell population abundance and this time-to-event variable. In practice, a frequent problem with time-to-even variables is that not all events are observed. Instead it is possible that only a minimum time is known, with other words those values are (right) censored. If standard analysis would be used (such as testDA_GLMM or testDA_edgeR from diffcyt) a bias is introduced in the parameter estimation. The simplest workaround is to exclude all incomplete samples, which is known under the name of complete case analysis or listwise deletion. But there are cases, such as high censoring rates or non random censoring, where this approach performs unfavorably, since it can greatly reduce statistical power or introduce a bias.
Normally, when a variable is censored, association testing can be done using classical survival analysis methods (e.g. Cox proportional hazard model). But in the context of differential abundance analysis in diffcyt the censored variable is not the response but a covariate (in contrary to classical survival analyis). The approach used by censcyt to overcome this hurdle is multiple imputation (Rubin 1988) which allows the use of the standard differential analysis methods at the cost of increased runtimes. Multiple imputation consists of the three steps Imputation, Analysis and Pooling which are shortly explained in the context of censcyt.
In the first step (Imputation) multiple complete data sets are generated by replacing the missing (or censored) values by a random draw from a predictive distribution (See Table 1). Then in the second step (Analysis) a GLMM is fitted on each imputed dataset where the cell population counts are modeled as the response and the censored (or rather imputed) variable (e.g. survival time) as a covariate. In the last step (Pooling) the results from the second step are combined using Rubins’s rules (Rubin 1988) which take into account the additional uncertainties from the imputation.
The three main available imputation methods in the function testDA_censoredGLMM are listed in Table 1.
| Method | Description | Abbr. |
|---|---|---|
| Complete case analysis | Incomplete samples are removed. | cc |
| Risk set imputation (Taylor et al. 2002) | Censored values are replace by a random value from its risk set. | rs |
| Kaplan-Meier imputation (Taylor et al. 2002) | Censored values are replaced by a random draw from a survival function that has been fit on the risk set of the respective censored value. | km |
| Mean residual life imputation (Conditional multiple imputation, Atem (2017)) | Bootstrapping of the incomplete data and replacement of the censored values by adding the mean residual life (the expected time until event given the censoring time). | mrl |
If the largest value is censored, the survival function does not reach its theoretical minimum of 0, which would leave some replacements undefined. Multiple options to handle this problem are implemented for use in Kaplan-Meier imputation (Table 2).
| Description | Abbr. |
|---|---|
| Set the largest value as observed. (default) | km |
| Model the tail of the survival function as an exponential distribution. (Moeschberger and Klein 1985) | km_exp |
| Model the tail of the survival function as a weibull distribution. (Moeschberger and Klein 1985) | km_wei |
| Replace all censored values larger than the largest observed value with expected values according to order statistics. (Moeschberger and Klein 1985) | km_os |
2 Data
First, load the necessary packages. Important to note is that loading censcyt will
also load diffcyt.
suppressPackageStartupMessages({
library(censcyt)
library(ggplot2)
library(SummarizedExperiment)
library(tidyr)
})The main input for DA analysis consists of a cluster times sample matrix of cell population abundances. How to obtain those counts is explained in the diffcyt vignette and will not be further considered here. Instead, for illustrative purposes, a dataset is simulated. In this case, the censored covariate is modeled according to an exponential distribution. To obtain censoring, an additional variable, the censoring time, is simulated as well and the actual measured value is the minimum of those two variables.
set.seed(1234)
nr_samples <- 50
nr_clusters <- 6
# the covariate is simulated from an exponential distribution:
X_true <- rexp(nr_samples)
# the censoring time is also sampled from an exponential distribution:
C <- rexp(nr_samples)
# the actual observed value is the minimum of the two:
X_obs <- pmin(X_true,C)
# additionally, we have the event indicator:
Event_ind <- ifelse(X_true>C,0,1)
# proportion censored:
(length(Event_ind)-sum(Event_ind))/length(Event_ind)## [1] 0.52Next, the cell counts are modeled according to a dirichlet-multinomial distribution where some clusters have an association with a covariate. The function simulate_multicluster can be used for this.
# number of cells per sample
sizes <- round(runif(nr_samples,1e3,1e4))
# alpha parameters of the dirichlet-multinomial distribution.
alphas <- runif(nr_clusters,1,10)
# main simulation function
simulation_output <- simulate_multicluster(alphas = alphas,
sizes = sizes,
covariate = X_obs,
nr_diff = 2,
return_summarized_experiment = TRUE,
slope = list(0.9))
# counts as a SummarizedExperiment object
d_counts <- simulation_output$countsTo get a sense of how this simulated data looks like we can plot the relative cell population abundance for each sample vs. the censored covariate per cell population.
# vector indicating if a cluster has a modeled association or not
is_diff_cluster <- ifelse(is.na(simulation_output$row_data$paired),FALSE,TRUE)
# first convert to proportions:
proportion <- as.data.frame(t(apply(assay(d_counts),2,function(x)x/sum(x))))
names(proportion) <- paste0("Nr: ",names(proportion), " - ",
ifelse(is_diff_cluster,"DA","non DA"))
# then to long format for plotting
counts_long <-
pivot_longer(proportion,
cols= seq_len(nr_clusters),
names_to = "cluster_id",
values_to = "proportion")
# add observed (partly censored) variable
counts_long$X_obs <- rep(X_obs,each=nr_clusters)
# add event indicator
counts_long$Event_ind <- factor(rep(Event_ind,each=nr_clusters),
levels = c(0,1),
labels=c("censored","observed"))
ggplot(counts_long) +
geom_point(aes(x=X_obs,y=proportion,color=Event_ind)) +
facet_wrap(~cluster_id)
Figure 1: Relative cell population abundance vs. survival time (X_obs)
Event_ind indicates if X_obs is censored or observed.
3 Set up meta-data
Next step is to set up the meta data, i.e. the covariates that are used in the testing. This is done the same way as in diffcyt but additionally the event indicator for the censored covariate has to be given, since censored variables are represented by two values, the measured value itself and an indicator if the measured value is observed (1) or censored (0).
# all covariates and sample ids
experiment_info <- data.frame(sample_id = seq_len(nr_samples),
X_obs = X_obs,
Event_ind = Event_ind)The createFormula function is similar to the one used in diffcyt but additionally the argument event_indicator has to be specified by giving the respective column name in experiment_info. If multiple covariates are given in argument cols_fixed then event_indicator is associated with the first element given to cols_fixed.
da_formula <- createFormula(experiment_info = experiment_info,
cols_fixed = "X_obs",
cols_random = "sample_id",
event_indicator = "Event_ind")
# also create contrast matrix for testing
contrast <- matrix(c(0,1))4 Differential testing
The main part is the differential testing which is done using the function testDA_censoredGLMM. It consists of the same arguments as the non-censored version testDA_GLMM (from diffcyt) with two additional arguments specific to multiple imputation. The first additional argument is mi_reps (multiple imputation repetitions) which is the number of multiple imputation steps performed. Unfortunately there is no clear rule as to how many imputations are needed and only rules of thumb are available (Buuren 2018). In general, more imputations are needed if the censoring rate is high and in applications where high statistical power is needed (Buuren 2018). One ‘rule’ is the quadratic rule of Hippel (2018) which can be shown if verbose is set to TRUE.
The second argument is imputation_method which is used to specify the imputation method. See Table 1 and Table 2.
# test with 50 repetitions with method risk set imputation (rs)
da_out <- testDA_censoredGLMM(d_counts = d_counts,
formula = da_formula,
contrast = contrast,
mi_reps = 30,
imputation_method = "km",
verbose = TRUE,
BPPARAM=BiocParallel::MulticoreParam(workers = 1))## 26 of 50 values are censored## suggested number of imputations: 54To look at the results we can use diffcyt’s topTable function:
topTable(da_out)## DataFrame with 6 rows and 3 columns
## cluster_id p_val p_adj
## <character> <numeric> <numeric>
## 3 3 0.000554004 0.00332403
## 6 6 0.004618430 0.01385529
## 1 1 0.780837186 0.99966962
## 2 2 0.970764737 0.99966962
## 4 4 0.999669622 0.99966962
## 5 5 0.679316197 0.99966962# compare to actual differential clusters:
which(!is.na(simulation_output$row_data$paired))## [1] 3 65 Wrapper function
Instead of using the single functions as shown above, it is possible use the wrapper
function censcyt which is the analog to the diffcyt wrapper. First, some expression data is simulated, i.e. a more realistic starting position than in the example above.
# Function to create expression data (one sample)
fcs_sim <- function(n = 2000, mean = 0, sd = 1, ncol = 10, cof = 5) {
d <- matrix(sinh(rnorm(n*ncol, mean, sd)) * cof,ncol=ncol)
# each marker is heavily expressed in one population, i.e. this should result
# in 'ncol' clusters
for(i in seq_len(ncol)){
d[seq(n/ncol*(i-1)+1,n/ncol*(i)),i] <-
sinh(rnorm(n/ncol, mean+2, sd)) * cof
}
colnames(d) <- paste0("marker", sprintf("%02d", 1:ncol))
d
}
# create multiple expression data samples with DA signal
comb_sim <- function(X, nr_samples = 10, mean = 0, sd = 1, ncol = 10, cof = 5){
# Create random data (without differential signal)
d_input <- lapply(seq_len(nr_samples),
function(i) fcs_sim(mean = mean,
sd = sd,
ncol = ncol,
cof = cof))
# Add differential abundance (DA) signal
for(i in seq_len(nr_samples)){
# number of cells in cluster 1
n_da <- round((plogis(X_true[i]))*300)
# set to random expression
d_input[[i]][seq_len(n_da), ] <-
matrix(sinh(rnorm(n_da*ncol, mean, sd)) * cof, ncol=ncol)
# increase expression for cluster 1
d_input[[i]][seq_len(n_da) ,1] <-
sinh(rnorm(n_da, mean + 2, sd)) * cof
}
d_input
}
# set parameter and simulat
d_input <- comb_sim(X_true, nr_samples = 50)The objects experiment_info, da_formula and contrast are the same as specified above, but additionally information about the measured markers has to be supplied.
marker_info <- data.frame(
channel_name = paste0("channel", sprintf("%03d", seq_len(10))),
marker_name = paste0("marker", sprintf("%02d", seq_len(10))),
marker_class = factor(c(rep("type", 10)),
levels = c("type", "state", "none")),
stringsAsFactors = FALSE
)Finally the wrapper function can be run. Argument mi_reps is to specify the number of repetitions in the multiple imputation and imputation_method allows to set the imputation method. To speed up the computation, a BiocParallelParam object (e.g. MulticoreParam, from package BiocParallel) can be given to the argument BPPARAM to parallelize parts of the computation.
# Run wrapper function
out_DA <- censcyt(d_input,
experiment_info,
marker_info,
formula = da_formula,
contrast = contrast,
meta_clustering = TRUE,
meta_k = 10,
seed_clustering = 123,
verbose = FALSE,
mi_reps = 10,
imputation_method = "km",
BPPARAM=BiocParallel::MulticoreParam(workers = 1))## FlowSOM clustering completed in 4.3 seconds# Display results for top DA clusters
topTable(out_DA, format_vals = TRUE)## DataFrame with 10 rows and 3 columns
## cluster_id p_val p_adj
## <character> <numeric> <numeric>
## 3 3 0.00422 0.0242
## 8 8 0.00483 0.0242
## 10 10 0.25100 0.8370
## 1 1 0.98000 0.9800
## 2 2 0.59100 0.9800
## 4 4 0.97500 0.9800
## 5 5 0.93800 0.9800
## 6 6 0.66000 0.9800
## 7 7 0.46400 0.9800
## 9 9 0.89800 0.9800# Plot heatmap for DA tests
plotHeatmap(out_DA,
analysis_type = "DA",
sample_order = order(X_true))
Session info
## R version 4.1.3 (2022-03-10)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.4 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] tidyr_1.2.0 SummarizedExperiment_1.24.0
## [3] Biobase_2.54.0 GenomicRanges_1.46.1
## [5] GenomeInfoDb_1.30.1 IRanges_2.28.0
## [7] S4Vectors_0.32.4 BiocGenerics_0.40.0
## [9] MatrixGenerics_1.6.0 matrixStats_0.61.0
## [11] ggplot2_3.3.5 censcyt_1.2.1
## [13] diffcyt_1.14.0 BiocStyle_2.22.0
##
## loaded via a namespace (and not attached):
## [1] backports_1.4.1 circlize_0.4.14
## [3] plyr_1.8.7 igraph_1.2.11
## [5] ConsensusClusterPlus_1.58.0 splines_4.1.3
## [7] flowCore_2.6.0 BiocParallel_1.28.3
## [9] listenv_0.8.0 scattermore_0.8
## [11] TH.data_1.1-0 digest_0.6.29
## [13] foreach_1.5.2 htmltools_0.5.2
## [15] magick_2.7.3 fansi_1.0.3
## [17] magrittr_2.0.2 CytoML_2.6.0
## [19] cluster_2.1.3 doParallel_1.0.17
## [21] aws.signature_0.6.0 limma_3.50.1
## [23] globals_0.14.0 ComplexHeatmap_2.10.0
## [25] RcppParallel_5.1.5 sandwich_3.0-1
## [27] flowWorkspace_4.6.0 cytolib_2.6.2
## [29] jpeg_0.1-9 colorspace_2.0-3
## [31] ggrepel_0.9.1 xfun_0.30
## [33] dplyr_1.0.8 crayon_1.5.1
## [35] RCurl_1.98-1.6 jsonlite_1.8.0
## [37] hexbin_1.28.2 graph_1.72.0
## [39] lme4_1.1-28 dirmult_0.1.3-5
## [41] survival_3.3-1 zoo_1.8-9
## [43] iterators_1.0.14 glue_1.6.2
## [45] polyclip_1.10-0 gtable_0.3.0
## [47] zlibbioc_1.40.0 XVector_0.34.0
## [49] GetoptLong_1.0.5 DelayedArray_0.20.0
## [51] ggcyto_1.22.0 car_3.0-12
## [53] Rgraphviz_2.38.0 shape_1.4.6
## [55] abind_1.4-5 scales_1.1.1
## [57] mvtnorm_1.1-3 pheatmap_1.0.12
## [59] DBI_1.1.2 edgeR_3.36.0
## [61] rstatix_0.7.0 Rcpp_1.0.8.3
## [63] clue_0.3-60 FlowSOM_2.2.0
## [65] httr_1.4.2 RColorBrewer_1.1-2
## [67] ellipsis_0.3.2 mice_3.14.0
## [69] pkgconfig_2.0.3 XML_3.99-0.9
## [71] farver_2.1.0 sass_0.4.1
## [73] locfit_1.5-9.5 utf8_1.2.2
## [75] labeling_0.4.2 reshape2_1.4.4
## [77] tidyselect_1.1.2 rlang_1.0.2
## [79] munsell_0.5.0 tools_4.1.3
## [81] cli_3.2.0 generics_0.1.2
## [83] broom_0.7.12 aws.s3_0.3.21
## [85] evaluate_0.15 stringr_1.4.0
## [87] fastmap_1.1.0 yaml_2.3.5
## [89] knitr_1.38 purrr_0.3.4
## [91] future_1.24.0 RBGL_1.70.0
## [93] nlme_3.1-157 xml2_1.3.3
## [95] compiler_4.1.3 curl_4.3.2
## [97] png_0.1-7 ggsignif_0.6.3
## [99] tibble_3.1.6 tweenr_1.0.2
## [101] bslib_0.3.1 stringi_1.7.6
## [103] highr_0.9 forcats_0.5.1
## [105] lattice_0.20-45 Matrix_1.4-1
## [107] nloptr_2.0.0 vctrs_0.3.8
## [109] furrr_0.2.3 pillar_1.7.0
## [111] lifecycle_1.0.1 BiocManager_1.30.16
## [113] jquerylib_0.1.4 GlobalOptions_0.1.2
## [115] data.table_1.14.2 bitops_1.0-7
## [117] colorRamps_2.3 R6_2.5.1
## [119] latticeExtra_0.6-29 bookdown_0.25
## [121] gridExtra_2.3 RProtoBufLib_2.6.0
## [123] parallelly_1.30.0 codetools_0.2-18
## [125] boot_1.3-28 MASS_7.3-56
## [127] assertthat_0.2.1 rjson_0.2.21
## [129] withr_2.5.0 broom.mixed_0.2.9.3
## [131] multcomp_1.4-18 GenomeInfoDbData_1.2.7
## [133] parallel_4.1.3 ncdfFlow_2.40.0
## [135] grid_4.1.3 ggpointdensity_0.1.0
## [137] minqa_1.2.4 rmarkdown_2.13
## [139] carData_3.0-5 Rtsne_0.15
## [141] ggpubr_0.4.0 ggnewscale_0.4.7
## [143] ggforce_0.3.3 base64enc_0.1-3References
Atem, Folefac D. 2017. “Linear Regression Model with a Randomly Censored Predictor: Estimation Procedures.” Biostatistics and Biometrics Open Access Journal 1 (2). https://doi.org/10.19080/bboaj.2017.01.555556.
Buuren, Stef van. 2018. Flexible Imputation of Missing Data, Second Edition. https://doi.org/10.1201/9780429492259.
Hippel, Paul T. von. 2018. “How Many Imputations Do You Need? A Two-stage Calculation Using a Quadratic Rule.” Sociological Methods & Research. https://doi.org/10.1177/0049124117747303.
Moeschberger, M. L., and John P. Klein. 1985. “A Comparison of Several Methods of Estimating the Survival Function When There is Extreme Right Censoring.” Biometrics 41 (1): 253. https://doi.org/10.2307/2530660.
Rubin, Donald B. 1988. “An Overview of Multiple Imputation.” Proceedings of the Survey Research Methods Section of the American Statistical Association.
Taylor, Jeremy M.G., Kristine L. Cooper, John T. Wei, Aruna V. Sarma, Trivellore E. Raghunathan, and Steve G. Heeringa. 2002. “Use of multiple imputation to correct for nonresponse bias in a survey of urologic symptoms among African-American men.” American Journal of Epidemiology. https://doi.org/10.1093/aje/kwf110.
Weber, Lukas M., Malgorzata Nowicka, Charlotte Soneson, and Mark D. Robinson. 2019. “diffcyt: Differential discovery in high-dimensional cytometry via high-resolution clustering.” Communications Biology. https://doi.org/10.1038/s42003-019-0415-5.