Hooke can identify differentially abundant cell types for a given contrast. The underlying PLN model mimics the standard (G)LM-like interface of R::stats. This allows users to fit a multivariate Poisson lognormal model after correcting for effects of offsets and covariates. This framework allows users to perform multivariate statistical regression to describe how perturbations alter the relative abundances of each cell state. For example, you could model the effects of genotype (WT vs MT) while controlling for experimental batch.
This data was published in Hasegawa, Franks, et al. bioRxiv. +The authors performed a longitudinal transcriptomic analysis using scRNA-sequencing of the lungs from pre- and post-silica challenged mice.
+The data includes 12 whole lung samples across pre- and post- intratracheal silica. 35 unique cell states were identified using highly and specifically expressed marker genes. For simplicity, we are splitting samples into 2 levels: exposed and not exposed. The goal is to identify which cell types are differentially abundant post-silica exposure.
+cds = readRDS("silicosis_cds.cds")
+
+# for simplicity, we are lumping together pre and post i.t. silica
+colData(cds)$exposed = ifelse(colData(cds)$Timepoint == 0, "not exposed", "exposed")
+colData(cds)$Rep = as.factor(colData(cds)$Rep)
+
+plot_cells(cds, color_cells_by = "fine_annotation", label_groups_by_cluster = F)
+
+
The first step of Hooke is to create a cell_count_set object. This requires three inputs:
cds - a monocle3 cell_data_set objectsample group - A column in colData(cds) that specifies how cells are grouped into samplescell group - A column in colData(cds) that specifies how cells are grouped into types or states (e.g. cluster)ccs = new_cell_count_set(cds,
+ sample_group = "ID",
+ cell_group = "fine_annotation")
+The next step is to fit a cell_count_model given a model formula. This function works with the standard R formula notations.
A cell count model requires the following inputs:
+ccs - a Hooke cell_count_set objectmain_model_formula_str - A character string specifying the model of cell abundances across samples, where terms refer to columns in colData(ccs). Put main effects here.nuisance_model_formula_str - A character string specifying the model of cell abundances across samples. Put nuisance effects here.In this case we are interested in comparing silica exposed mouse lungs to unexposed lungs.
+ccm = new_cell_count_model(ccs,
+ main_model_formula_str = "~ exposed")
+Next we can predict fold change values for a given contrast of interest using estimate_abundances().
ccm - a Hooke cell_count_model objectnewdata - A tibble of variables used for the prediction.min_log_abund - numeric Minimum log abundance value.cond_exp = estimate_abundances(ccm, tibble::tibble(exposed = "exposed"))
+cond_not_exp = estimate_abundances(ccm, tibble::tibble(exposed = "not exposed"))
+
+cond_not_exp %>% head()
+| exposed | +cell_group | +log_abund | +log_abund_se | +log_abund_sd | +
|---|---|---|---|---|
| not exposed | +Endothelial Arterial | +2.87 | +0.196 | +0.559 | +
| not exposed | +Endothelial Capillary | +5.31 | +0.0591 | +0.559 | +
| not exposed | +Endothelial Venous | +3.09 | +0.179 | +0.559 | +
| not exposed | +Lymphatic Endothelial | +2.63 | +0.213 | +0.559 | +
| not exposed | +Aerocytes | +4.06 | +0.110 | +0.559 | +
| not exposed | +Fibroblasts | +5.26 | +0.0607 | +0.559 | +
And then compare two estimates of cell abundances from a Hooke model to identify significant changes with compare_abundances(). In this case, exposed vs not exposed.
ccm - a Hooke cell_count_model objectcond_x - A cell type abundance estimate from estimate_abundances().cond_y - A cell type abundance estimate from estimate_abundances().cond_ne_v_e_tbl = compare_abundances(ccm, cond_not_exp, cond_exp)
+
+cond_ne_v_e_tbl %>% select(cell_group, exposed_x, exposed_y,
+ delta_log_abund, delta_log_abund_se, delta_q_value)
+
+| cell_group | +delta_log_abund | +delta_log_abund_se | +delta_q_value | +
|---|---|---|---|
| Endothelial Arterial | +-0.520 | +0.268 | +0.225 | +
| Endothelial Capillary | +-0.400 | +0.0799 | +0.0109 | +
| Endothelial Venous | +-0.332 | +0.237 | +0.356 | +
| Lymphatic Endothelial | +0.0650 | +0.283 | +0.925 | +
| Aerocytes | +-0.351 | +0.148 | +0.143 | +
| Fibroblasts | +-0.477 | +0.0877 | +0.0109 | +
Finally, we can color our UMAP by the fold changes we estimate in compare_abundances() using the plot_contrast() function. We can select fold changes that are signifiant given a q_value_thresh.
ccm - A Hooke cell_count_model object.cond_b_vs_a_tbl - A data frame from compare_abundances().log_abundance_thresh - numeric Select cell groups by log abundance.q_value_thresh - Remove contrasts whose change in q-value exceeds q_value_thresh.plot_contrast(ccm, cond_ne_v_e_tbl, q_value_thresh = 0.05)
+
If your data contains multiple experimental batch, Hooke supports batch correction with the use of a nuisance model string term.
+ccm_rep = new_cell_count_model(ccs,
+ main_model_formula_str = "~ exposed",
+ nuisance_model_formula_str = "~ Rep")
+
+cond_exp_rep = estimate_abundances(ccm_rep, tibble::tibble(exposed = "exposed", Rep = "1"))
+cond_not_exp_rep = estimate_abundances(ccm_rep, tibble::tibble(exposed = "not exposed", Rep = "1"))
+cond_ne_v_e_tbl_rep = compare_abundances(ccm_rep, cond_not_exp_rep, cond_exp_rep)
+
+plot_contrast(ccm_rep, cond_ne_v_e_tbl_rep, q_value_thresh = 0.05)
+Transitional interstitial macrophages and regulatory T cells no longer come up as significant if we control for replicate.
+
Hooke can be used to analyze data with multiple time points across perturbations, such as the data found in +Saunders, Srivatsan, et al. Nature, in press (2023). This study in includes ~3 million single cells across almost 2000 individual barcoded embryos. It includes 19 timepoints (18-96 hpf) and 23 genetic loss of function experiments. For more information see the ZSCAPE website.
+This vignette focuses on the cranial sensory ganglia subset. Cranial sensory neurons transmit information from the head, ear, heart and viscera. Around 30,000 cranial sensory neurons (~20 cells/embryo) were identified. These cells were annotated by comparing branch specific gene expression with published expression data.
+
Foxi1 and phox2a are two transcription factors important for cranial sensory ganglia development. Here we subset the cds to focus on those perturbations. We also include all the controls in our cell count set construction.
+control_ids = c("ctrl-uninj", "ctrl-inj", "ctrl-hgfa", "ctrl-met", "ctrl-mafba", "ctrl-noto", "ctrl-tbx16")
+cds = cds[, colData(cds)$gene_target %in% c("phox2a", "foxi1", control_ids)]
+
+colData(cds)$perturbation = ifelse(colData(cds)$gene_target %in% control_ids,
+ "control",
+ colData(cds)$gene_target)
+
+ccs = new_cell_count_set(cds,
+ sample_group = "embryo",
+ cell_group = "cell_type_sub")
+The functionbuild_interval_formula() builds a model formula for time series models based on the range of the data. This is a utility function that puts the knots in reasonable positions based on the range of the data. It takes the following as input:
ccs- a Hooke cell_count_set object. num_breaks - specifies degrees of freedominterval_var - name of the column that the interval should be based on interval_start- defaults to the min(interval_var) valueinterval_stop - defaults to the max(interval_var) valueNote: if using this function or specifying your own spline, you must load library(spline). Putting splines::ns() in your model formula will result in an error as the model reads it as an interaction term.
start_time = 18
+stop_time = 72
+time_formula = build_interval_formula(ccs, num_breaks = 3, interval_start = 18, interval_stop = 72)
+The resulting time formula is ~ ns(timepoint , knots= c(45))
We can fit a cell_count_model using this time formula and a perturbation term...
ccm = new_cell_count_model(ccs,
+ main_model_formula_str = paste0("perturbation +", time_formula))
+
+# predict for 48 hpf
+cond_wt = estimate_abundances(ccm, tibble(timepoint = 48, perturbation = "control"))
+cond_phox2a = estimate_abundances(ccm, tibble(timepoint = 48, perturbation = "phox2a"))
+cond_foxi1 = estimate_abundances(ccm, tibble(timepoint = 48, perturbation = "foxi1"))
+
+wt_v_phox2a_tbl = compare_abundances(ccm, cond_wt, cond_phox2a)
+wt_v_foxi1_tbl = compare_abundances(ccm, cond_wt, cond_foxi1)
+... and plot the abundance changes.
+plot_contrast(ccm, wt_v_phox2a_tbl, x=1, y=3, q_value_thresh = 0.05)
+
We can use Hooke to view when cell types are at their peak abundance in the dataset. Here we are just looking at wild-type data and estimating abundances over a time interval.
+
+wt_cds = cds[, colData(cds)$gene_target %in% control_ids]
+wt_ccs = new_cell_count_set(wt_cds,
+ sample_group = "embryo",
+ cell_group = "cell_type_sub")
+
+colData(wt_ccs)$timepoint = as.numeric(colData(wt_ccs)$timepoint)
+wt_ccm = new_cell_count_model(wt_ccs,
+ main_model_formula_str = "ns(timepoint, df=3)")
+
+wt_timepoint_pred_df = estimate_abundances_over_interval(wt_ccm,
+ interval_start=18,
+ interval_stop=72,
+ interval_col="timepoint",
+ interval_step=2)
+
+log_abund_detection_thresh=-3
+
+ggplot(wt_timepoint_pred_df, aes(x = timepoint)) +
+ geom_line(aes(y = exp(log_abund) + exp(log_abund_detection_thresh))) +
+ facet_wrap(~cell_group, scales="free_y") + monocle3:::monocle_theme_opts()
+
+
If the data was collected in multiple batches, you can also include a nuisance_model_formula_str and plot each experimental batch's kinetics.
+wt_expt_ccm = new_cell_count_model(wt_ccs,
+ main_model_formula_str = "ns(timepoint, df=3)",
+ nuisance_model_formula_str = "~ expt")
+
+batches = data.frame(batch = unique(colData(wt_ccs)$expt))
+batches = batches %>% mutate(tp_preds = purrr::map(.f = function(batch) {
+ estimate_abundances_over_interval(wt_expt_ccm,
+ start_time,
+ stop_time,
+ knockout=FALSE,
+ interval_col="timepoint",
+ interval_step=2,
+ expt = batch)
+}, .x=batch))
+
+wt_timepoint_pred_df = batches %>% select(tp_preds) %>% tidyr::unnest(tp_preds)
+
+ggplot(wt_timepoint_pred_df, aes(x = timepoint)) +
+ geom_line(aes(y = exp(log_abund) + exp(log_abund_detection_thresh), color=expt)) +
+ facet_wrap(~cell_group, scales="free_y", nrow = 2) + monocle3:::monocle_theme_opts() +
+ ggtitle("wild-type kinetics by expt")
+
+
We can also plot a given perturbation's kinetics against the wild type.
+
+foxi1_cds = cds[, colData(cds)$gene_target %in% c("foxi1", control_ids)]
+foxi1_ccs = new_cell_count_set(foxi1_cds,
+ sample_group = "embryo",
+ cell_group = "cell_type_sub")
+
+colData(foxi1_ccs)$timepoint = as.numeric(colData(foxi1_ccs)$timepoint)
+
+
+start_time = 18
+stop_time = 72
+
+time_formula = build_interval_formula(foxi1_ccs, num_breaks = 3, interval_start = 18, interval_stop = 72)
+
+foxi1_ccm = new_cell_count_model(foxi1_ccs,
+ main_model_formula_str = paste0("perturbation +", time_formula))
+
+wt_timepoint_pred_df = estimate_abundances_over_interval(foxi1_ccm,
+ interval_start=start_time,
+ interval_stop=stop_time,
+ interval_col="timepoint",
+ interval_step=2,
+ perturbation = "control")
+
+ko_timepoint_pred_df = estimate_abundances_over_interval(foxi1_ccm,
+ interval_start=start_time,
+ interval_stop=stop_time,
+ interval_col="timepoint",
+ interval_step=2,
+ perturbation = "foxi1")
+timepoints = seq(start_time, stop_time, 2)
+perturb_vs_wt_nodes = tibble(t1=timepoints) %>%
+ mutate(comp_abund = purrr::map(.f = compare_ko_to_wt_at_timepoint,
+ .x = t1,
+ perturbation_ccm=foxi1_ccm,
+ interval_col="timepoint",
+ wt_pred_df = wt_timepoint_pred_df,
+ ko_pred_df = ko_timepoint_pred_df)) %>%
+ tidyr::unnest(comp_abund)
+
+
+ggplot(perturb_vs_wt_nodes, aes(x = t1)) +
+ geom_line(aes(y = exp(log_abund_x) + exp(log_abund_detection_thresh), linetype = "Wild-type")) +
+ geom_line(aes(y = exp(log_abund_y) + exp(log_abund_detection_thresh), linetype = "Knockout")) +
+ ggh4x::stat_difference(aes(ymin = exp(log_abund_x)+exp(log_abund_detection_thresh), ymax = exp(log_abund_y) +exp(log_abund_detection_thresh)), alpha=0.3) +
+ facet_wrap(~cell_group, scales="free_y", nrow=2) + monocle3:::monocle_theme_opts() + ggtitle("foxi1 kinetics")
+
+
+