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 @ref(tab:methods)). 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 @ref(tab:methods).

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 @ref(tab:tailmethods)).

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

Next, 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$counts

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

Relative cell population abundance vs. survival time (X_obs). Event_ind indicates if X_obs is censored or observed.

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

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 @ref(tab:methods) and Table @ref(tab:tailmethods).

# 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:  54

To 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.00148338 0.00890028
## 6           6 0.00987323 0.02961970
## 1           1 0.57299042 0.89086778
## 5           5 0.59391185 0.89086778
## 4           4 0.77492552 0.92991063
## 2           2 0.98172528 0.98172528

# compare to actual differential clusters:
which(!is.na(simulation_output$row_data$paired))

## [1] 3 6

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 3.5 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>
## 8            8  0.000377   0.00377
## 3            3  0.001770   0.00883
## 1            1  0.833000   0.96600
## 2            2  0.557000   0.96600
## 4            4  0.869000   0.96600
## 5            5  0.616000   0.96600
## 6            6  0.835000   0.96600
## 7            7  0.535000   0.96600
## 10          10  0.313000   0.96600
## 9            9  0.973000   0.97300

# Plot heatmap for DA tests
plotHeatmap(out_DA, 
            analysis_type = "DA", 
            sample_order = order(X_true))