3. Case Study: A Disease-Associated Monocyte State in COVID-19
23 September 2026
Package
scDiagnostics 1.7.15
1 Purpose
This vignette walks through the project-detect-characterize workflow introduced in vignette 1 on a real disease case study: PBMC scRNA-seq data from healthy donors (reference) and donors with severe COVID-19 (query), from Stephenson et al. (2021). The goal is to find out whether CD14 monocytes in the severe-COVID query look like a distinct state relative to the healthy reference, and if so, what distinguishes them.
library(scDiagnostics)
library(SingleCellExperiment)
set.seed(100)2 The data
covid_reference_data (healthy donors) and covid_query_data (severe
COVID-19 donors) are downsampled subsets of the Stephenson et al. (2021)
PBMC atlas, restricted to 5 shared cell types and a gene panel that
always includes a 25-gene interferon-response signature (Yoshida et
al.); see ?covid_reference_data for the full processing details.
data("covid_reference_data")
data("covid_query_data")
table(covid_reference_data$author_cell_type_merged)
#>
#> B cell CD14 mono CD4 T CD8 T NK_16hi
#> 180 180 180 180 180
table(covid_query_data$azimuth_celltype_l1_merged)
#>
#> B cell CD14 mono CD4 T CD8 T NK_16hi
#> 220 450 220 220 220The reference’s cell type column (author_cell_type_merged) reflects
the original authors’ annotation; the query’s (azimuth_celltype_l1_merged)
comes from Azimuth reference mapping. Both were computed independently
of scDiagnostics - we are auditing an annotation transfer that has
already happened, not producing one.
3 Step 1: project
plotCellTypePCA() projects the query onto the reference’s PCA space
and compares the distributions of each cell type along the leading PCs:
shared_cell_types <- c("CD14 mono", "CD4 T", "CD8 T", "B cell", "NK_16hi")
plotCellTypePCA(
query_data = covid_query_data,
reference_data = covid_reference_data,
cell_types = shared_cell_types,
query_cell_type_col = "azimuth_celltype_l1_merged",
ref_cell_type_col = "author_cell_type_merged",
pc_subset = 1:3)
#> Picking joint bandwidth of 0.304
#> Picking joint bandwidth of 0.589
#> Picking joint bandwidth of 0.491
4 Step 2: detect
Focusing specifically on CD14 monocytes, detectAnomaly() builds an
Isolation Forest on the reference’s PCA projection and scores how
anomalous each query cell looks relative to it:
anomaly_output <- detectAnomaly(
reference_data = covid_reference_data,
query_data = covid_query_data,
ref_cell_type_col = "author_cell_type_merged",
query_cell_type_col = "azimuth_celltype_l1_merged",
cell_types = "CD14 mono",
pc_subset = 1:5,
n_tree = 500)
mean(anomaly_output[["CD14 mono"]]$query_anomaly)
#> [1] 0.2733333plot(anomaly_output, cell_type = "CD14 mono", pc_subset = 1:3, data_type = "query")
In this downsampled dataset, a substantial fraction of the query’s CD14 monocytes are flagged as anomalous relative to the healthy reference. Unlike the rare/withheld-cell-type scenarios in vignette 2, this is not necessarily a small, rare subpopulation - a disease process can plausibly shift a large fraction of a cell type’s expression profile, and that is a hypothesis worth checking directly rather than assuming anomaly detection here means the same thing it did there.
5 Step 3: characterize
calculateGeneShifts() tests each gene in a specified panel for a
distributional shift between reference and query, optionally comparing
only the non-anomalous reference cells against the anomalous query
cells (anomaly_comparison = TRUE). We focus this on the 25-gene
Yoshida et al. interferon-response signature already included in the
gene panel of both objects:
yoshida_ifn_signature <- c(
"BST2", "CMPK2", "EIF2AK2", "EPSTI1", "HERC5", "IFI35", "IFI44L",
"IFI6", "IFIT3", "ISG15", "LY6E", "MX1", "MX2", "OAS1", "OAS2",
"PARP9", "PLSCR1", "SAMD9", "SAMD9L", "SP110", "STAT1", "TRIM22",
"UBE2L6", "XAF1", "IRF7")
gene_shifts <- calculateGeneShifts(
query_data = covid_query_data[yoshida_ifn_signature, ],
reference_data = covid_reference_data[yoshida_ifn_signature, ],
query_cell_type_col = "azimuth_celltype_l1_merged",
ref_cell_type_col = "author_cell_type_merged",
cell_types = "CD14 mono",
pc_subset = 1:5,
n_top_loadings = 25,
detect_anomalies = TRUE,
anomaly_comparison = TRUE)
#> Warning in check_numbers(x, k = k, nu = nu, nv = nv): more singular
#> values/vectors requested than available
head(gene_shifts$PC1[order(gene_shifts$PC1$p_adjusted), ], 10)
#> gene loading cell_type p_value mean_query mean_reference p_adjusted
#> 1 LY6E 0.017244450 CD14 mono 0 2.7398890 0.80911272 0
#> 2 IFI6 -0.014406456 CD14 mono 0 2.6754521 0.47198159 0
#> 3 BST2 -0.012502318 CD14 mono 0 1.6340352 0.82511624 0
#> 4 IRF7 -0.011014061 CD14 mono 0 1.3113019 0.39981541 0
#> 5 MX2 -0.010216540 CD14 mono 0 0.9759931 0.28497303 0
#> 6 MX1 -0.008912552 CD14 mono 0 1.3944989 0.24070821 0
#> 7 IFI44L -0.007559699 CD14 mono 0 1.1655427 0.13011109 0
#> 8 ISG15 -0.006431524 CD14 mono 0 2.0172213 0.29153589 0
#> 9 EPSTI1 -0.004185579 CD14 mono 0 1.1093122 0.13667800 0
#> 10 IFIT3 -0.003400604 CD14 mono 0 0.6571842 0.04993755 0
#> significant
#> 1 TRUE
#> 2 TRUE
#> 3 TRUE
#> 4 TRUE
#> 5 TRUE
#> 6 TRUE
#> 7 TRUE
#> 8 TRUE
#> 9 TRUE
#> 10 TRUEThe genes with the smallest adjusted p-values here (e.g. IFI6, LY6E,
BST2) are all canonical interferon-stimulated genes, each with
substantially higher mean expression in the anomalous query cells than
in the non-anomalous reference cells. plot() visualizes this as a
heatmap of per-gene z-scores:
plot(gene_shifts, cell_type = "CD14 mono", pc_subset = 1:5,
plot_type = "heatmap", plot_by = "p_adjusted", n_genes = 15,
show_anomalies = TRUE, pseudo_bulk = TRUE, cluster_cols = TRUE)
or as fold-changes relative to the reference:
plot(gene_shifts, cell_type = "CD14 mono", pc_subset = 1:5,
plot_type = "barplot", plot_by = "p_adjusted", n_genes = 15,
show_anomalies = TRUE, pseudo_bulk = TRUE)
Together, these three steps give a concrete, checkable answer: yes,
CD14 monocytes in the severe-COVID query look different from the
healthy reference in PCA space, detectAnomaly() flags a large fraction
of them accordingly, and the genes distinguishing the flagged cells are
specifically interferon-response genes - consistent with a known
biological interferon-activated monocyte state in severe COVID-19,
rather than an artifact of the annotation transfer itself. The original
manuscript further shows this same interferon signature recovered
regardless of which of four independent annotation tools (Azimuth,
SingleR, CellTypist, scVI) produced the query labels; that
cross-tool comparison is not reproduced in this vignette.
6 R Session Info
R version 4.6.1 (2026-06-24)
Platform: x86_64-pc-linux-gnu
Running under: Ubuntu 24.04.5 LTS
Matrix products: default
BLAS: /home/biocbuild/bbs-3.24-bioc/R/lib/libRblas.so
LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.12.0 LAPACK version 3.12.0
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
time zone: America/New_York
tzcode source: system (glibc)
attached base packages:
[1] stats4 stats graphics grDevices utils datasets methods
[8] base
other attached packages:
[1] SingleCellExperiment_1.35.2 SummarizedExperiment_1.43.0
[3] Biobase_2.73.2 GenomicRanges_1.65.4
[5] Seqinfo_1.3.2 IRanges_2.47.5
[7] S4Vectors_0.51.10 BiocGenerics_0.59.12
[9] generics_0.1.4 MatrixGenerics_1.25.0
[11] matrixStats_1.5.0 scDiagnostics_1.7.15
[13] BiocStyle_2.41.0
loaded via a namespace (and not attached):
[1] gridExtra_2.3.1 rlang_1.3.0 magrittr_2.0.5
[4] clue_0.3-68 GetoptLong_1.1.1 scater_1.41.2
[7] otel_0.2.0 ggridges_0.5.7 compiler_4.6.1
[10] png_0.1-9 systemfonts_1.3.2 vctrs_0.7.3
[13] shape_1.4.6.1 crayon_1.5.3 pkgconfig_2.0.3
[16] fastmap_1.2.0 magick_2.9.1 XVector_0.53.0
[19] scuttle_1.23.2 labeling_0.4.3 rmarkdown_2.32
[22] ggbeeswarm_0.7.3 ragg_1.5.2 tinytex_0.61
[25] purrr_1.2.2 xfun_0.61 bluster_1.23.1
[28] cachem_1.1.0 beachmat_2.29.3 jsonlite_2.0.0
[31] DelayedArray_0.39.7 BiocParallel_1.47.0 irlba_2.3.7
[34] parallel_4.6.1 cluster_2.1.8.3 R6_2.6.1
[37] bslib_0.12.0 RColorBrewer_1.1-3 limma_3.99.0
[40] GGally_2.4.0 jquerylib_0.1.4 iterators_1.0.14
[43] Rcpp_1.1.2 bookdown_0.48 knitr_1.52
[46] Matrix_1.7-6 igraph_2.3.3 tidyselect_1.2.1
[49] dichromat_2.0-1 abind_1.4-8 yaml_2.3.12
[52] viridis_0.6.5 doParallel_1.0.17 codetools_0.2-20
[55] lattice_0.23-1 tibble_3.3.1 withr_3.0.3
[58] S7_0.2.2 evaluate_1.0.5 ggstats_0.14.0
[61] circlize_0.4.18 pillar_1.11.1 BiocManager_1.30.27
[64] foreach_1.5.2 ggplot2_4.0.3 scales_1.4.0
[67] RhpcBLASctl_0.23-42 glue_1.8.1 metapod_1.21.0
[70] tools_4.6.1 BiocNeighbors_2.7.3 ScaledMatrix_1.21.0
[73] locfit_1.5-9.12 scran_1.41.1 Cairo_1.7-0
[76] grid_4.6.1 tidyr_1.3.2 colorspace_2.1-3
[79] edgeR_4.99.6 beeswarm_0.4.0 BiocSingular_1.29.1
[82] vipor_0.4.7 cli_3.6.6 rsvd_1.0.5
[85] textshaping_1.0.5 S4Arrays_1.13.1 viridisLite_0.4.3
[88] ComplexHeatmap_2.29.0 dplyr_1.2.1 gtable_0.3.6
[91] isotree_0.6.1-5 sass_0.4.10 digest_0.6.39
[94] SparseArray_1.13.3 ggrepel_0.9.8 dqrng_0.4.1
[97] rjson_0.2.23 farver_2.1.2 htmltools_0.5.9
[100] lifecycle_1.0.5 GlobalOptions_0.1.4 statmod_1.5.2