suppressPackageStartupMessages({
library(fishpond)
library(SummarizedExperiment)
library(SingleCellExperiment)
library(scuttle)
library(scater)
library(scran)
library(dplyr)
library(tidyr)
library(ggplot2)
library(msigdbr)
library(cowplot)
library(ComplexHeatmap)
library(circlize)
library(scDblFinder)
library(Seurat)
library(sketchR)
library(tricycle)
library(org.Hs.eg.db)
library(parallel)
library(bluster)
library(BiocParallel)
library(igraph)
library(DT)
library(swissknife)
library(SingleR)
})Organoid scRNAseq timecourse
Load packages
Define parameters and reference files
## Sample annotation
sampleannot <- "scrnaseq_timecourse_samples.txt"
## Gene annotation (derived from GENCODE v38 gtf file)
t2g <- "reference/gencode.v38.annotation.tx2gene_full.rds"
## Cell cycle marker genes
## https://github.com/ventolab/HECA-Human-Endometrial-Cell-Atlas/tree/main/utils
s_genes <- "data/reproductivecellatlas/2023/S_genes_HECA.tsv"
g2m_genes <- "data/reproductivecellatlas/2023/G2M_genes_HECA.tsv"
## Tissue atlas
tissue_atlas_sce <- "../scRNAseq_invivo_endometrium_epithelial_cell_atlas/invivo_epithelial_cell_atlas.rds"Assemble gene annotation
geneannot <- readRDS(t2g) |>
dplyr::select(gene_id, seqnames, gene_type, gene_name) |>
dplyr::rename(chromosome = seqnames) |>
dplyr::distinct()
## Annotate ribosomal protein genes
## GO term "structural constituent of ribosome" (GO:0003735)
golist <- as.list(org.Hs.egGO2ALLEGS)
ribo_entrez <- unique(golist[["GO:0003735"]])
length(ribo_entrez)[1] 175tmp <- select(org.Hs.eg.db, keys = ribo_entrez, columns = "ENSEMBL",
keytype = "ENTREZID")'select()' returned 1:many mapping between keys and columnsribo_ensembl <- unique(na.omit(tmp[, "ENSEMBL"]))
length(ribo_ensembl)[1] 195summary(ribo_ensembl %in% sub("\\.[0-9]+$", "", geneannot$gene_id)) Mode FALSE TRUE
logical 33 162 geneannot$is_ribosomal_protein <-
sub("\\.[0-9]+$", "", geneannot$gene_id) %in% ribo_ensembl
geneannot$gene_type_ribo <- ifelse(
geneannot$gene_type == "protein_coding",
ifelse(geneannot$is_ribosomal_protein, "protein_coding_ribosome",
geneannot$gene_type),
geneannot$gene_type)Read data
## Get sample IDs
tbl <- read.delim(sampleannot)
DT::datatable(tbl, options = list(scrollX = TRUE, pageLength = 50))## Read alevin-fry data and combine
dat <- do.call(cbind, lapply(seq_len(nrow(tbl)), function(i) {
x <- fishpond::loadFry(
fryDir = file.path("data", tbl$dir[i], tbl$sample[i]),
outputFormat = list(counts = c("S", "A"), spliced = c("S", "A"),
unspliced = c("U"), total = c("U", "S", "A")))
stopifnot(all(!duplicated(colnames(x))))
rowData(x)$gene_name <-
geneannot$gene_name[match(rowData(x)$gene_ids, geneannot$gene_id)]
rowData(x)$chromosome <-
geneannot$chromosome[match(rowData(x)$gene_ids, geneannot$gene_id)]
rowData(x)$gene_type <-
geneannot$gene_type[match(rowData(x)$gene_ids, geneannot$gene_id)]
rowData(x)$gene_type_ribo <-
geneannot$gene_type_ribo[match(rowData(x)$gene_ids, geneannot$gene_id)]
## Keep only cells with at least 1,000 detected genes
x <- x[, colSums(assay(x, "total") > 0) > 1000]
## Add sample ID to colnames
x$sample <- tbl$sample[i]
colnames(x) <- paste0(tbl$sample[i], "_", colnames(x))
x
}))
datclass: SingleCellExperiment
dim: 60649 114374
metadata(0):
assays(4): counts spliced unspliced total
rownames(60649): ENSG00000223972.5 ENSG00000243485.5 ...
ENSG00000210194.1 ENSG00000210196.2
rowData names(5): gene_ids gene_name chromosome gene_type
gene_type_ribo
colnames(114374): 3684F1_ATGGTTGCAATTTCTC 3684F1_AAGTACCTCTCTGGTC ...
3773F9_ACAGCCGGTGAGTAAT 3773F9_TGACCCTCAACACGTT
colData names(2): barcodes sample
reducedDimNames(0):
mainExpName: NULL
altExpNames(0):## Set rownames to gene names where possible
rownames(dat) <- scater::uniquifyFeatureNames(
ID = rownames(dat),
names = rowData(dat)$gene_name
)Add sample annotations
## Add 'cell' column
dat$cell <- colnames(dat)
## Any other annotations
if (nrow(tbl) > 0 && (length(tbl$sample) == length(unique(tbl$sample)))) {
for (cn in colnames(tbl)) {
if (!all(is.na(tbl[[cn]])) && !(cn %in% colnames(colData(dat)))) {
colData(dat)[[cn]] <- tbl[[cn]][match(dat$sample, tbl$sample)]
}
}
}
## Create donorCol
dat$donorCol <- dat$donor
## Order timepoints
dat$timepoint <- c(d8 = "Horm Withdr 48h", d9_12h = "Post-breakdown 12h",
d9_24h = "Post-breakdown 24h", d10 = "Post-breakdown 48h",
d11 = "E2-diff 24h", d12 = "E2-diff 48h")[dat$timepoint]
dat$timepoint <- factor(dat$timepoint,
levels = c("Horm Withdr 48h", "Post-breakdown 12h",
"Post-breakdown 24h", "Post-breakdown 48h",
"E2-diff 24h", "E2-diff 48h"))
table(dat$timepoint, useNA = "ifany")
Horm Withdr 48h Post-breakdown 12h Post-breakdown 24h Post-breakdown 48h
24451 19557 17708 15825
E2-diff 24h E2-diff 48h
18038 18795 Check total count for each gene type and removed unwanted gene types
keepGeneTypes <- c("protein_coding", "IG_V_gene", "IG_C_gene", "IG_J_gene",
"TR_C_gene", "TR_J_gene", "TR_V_gene", "TR_D_gene",
"IG_D_gene", "Mt_tRNA", "Mt_rRNA")
pd <- data.frame(gene = rownames(dat),
gene_type_ribo = rowData(dat)$gene_type_ribo,
total_count = rowSums(counts(dat))) |>
group_by(gene_type_ribo) |>
summarize(nbr_genes = length(gene),
total_count = sum(total_count)) |>
mutate(keep = gene_type_ribo %in% keepGeneTypes) |>
arrange(desc(nbr_genes))
pd$gene_type_ribo <- factor(pd$gene_type_ribo, levels = unique(pd$gene_type_ribo))
p1 <- ggplot(pd, aes(x = gene_type_ribo, y = nbr_genes, fill = keep)) +
geom_col() +
labs(x = "", y = "Number of genes", fill = "keep") +
theme_cowplot() +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
p2 <- ggplot(pd, aes(x = gene_type_ribo, y = total_count, fill = keep)) +
geom_col() +
labs(x = "", y = "Total UMI count", fill = "keep") +
theme_cowplot() +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
plot_grid(p1, p2, ncol = 1)
dat <- dat[rowData(dat)$gene_type_ribo %in% keepGeneTypes, ]
dim(dat)[1] 20231 114374Add QC metrics
(MT <- grep("^MT-", rownames(dat), value = TRUE)) [1] "MT-TF" "MT-RNR1" "MT-TV" "MT-RNR2" "MT-TL1" "MT-ND1" "MT-TI"
[8] "MT-TM" "MT-ND2" "MT-TW" "MT-CO1" "MT-TD" "MT-CO2" "MT-TK"
[15] "MT-ATP8" "MT-ATP6" "MT-CO3" "MT-TG" "MT-ND3" "MT-TR" "MT-ND4L"
[22] "MT-ND4" "MT-TH" "MT-TS2" "MT-TL2" "MT-ND5" "MT-CYB" "MT-TT"
[29] "MT-TQ" "MT-TA" "MT-TN" "MT-TC" "MT-TY" "MT-TS1" "MT-ND6"
[36] "MT-TE" "MT-TP" dat <- scuttle::addPerCellQC(dat, assay.type = "counts", subsets = list(MT = MT))
dat <- scuttle::addPerFeatureQC(dat, assay.type = "counts")
dat$intronic_fraction <-
colSums(assay(dat, "unspliced")) / colSums(assay(dat, "total"))
rowData(dat)$intronic_fraction <-
rowSums(assay(dat, "unspliced")) / rowSums(assay(dat, "total"))QC plots
total_count_threshold <- 1000
detected_threshold <- 500
mt_percent_threshold <- 25
dat$keep <- (dat$total > total_count_threshold) &
(dat$detected > detected_threshold) &
(dat$subsets_MT_percent < mt_percent_threshold)
table(dat$keep)
FALSE TRUE
26943 87431 cdplot <- as.data.frame(colData(dat))
sid <- "sample"
stid <- "timepoint"
ggplot(cdplot, aes(x = .data[[sid]], y = sum, fill = .data[[stid]])) +
geom_hline(yintercept = total_count_threshold) +
geom_violin(scale = "width") + theme_minimal() +
labs(y = "Total UMI count per cell") +
theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5))
ggplot(cdplot, aes(x = .data[[sid]], y = sum, fill = .data[[stid]])) +
geom_hline(yintercept = total_count_threshold) +
geom_violin(scale = "width") + theme_minimal() +
scale_y_log10() +
labs(y = "Total UMI count per cell") +
theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5))
ggplot(cdplot, aes(x = .data[[sid]], y = detected, fill = .data[[stid]])) +
geom_hline(yintercept = detected_threshold) +
geom_violin(scale = "width") + theme_minimal() +
labs(y = "Number of detected genes per cell") +
theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5))
ggplot(cdplot, aes(x = .data[[sid]], y = subsets_MT_percent, fill = .data[[stid]])) +
geom_hline(yintercept = mt_percent_threshold) +
geom_violin(scale = "width") + theme_minimal() +
labs(y = "Percentage of MT count per cell") +
theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5))
ggplot(cdplot, aes(x = detected, y = subsets_MT_percent, color = keep)) +
geom_point(size = 0.5, alpha = 0.5) + theme_minimal() +
labs(x = "Number of detected genes per cell",
y = "Percentage of MT count per cell")
ggplot(cdplot, aes(x = detected, y = intronic_fraction, color = keep)) +
geom_point(size = 0.5, alpha = 0.5) + theme_minimal() +
labs(x = "Number of detected genes per cell",
y = "Intronic fraction per cell")
Filter and remove samples with too few cells
(tbl <- table(sample = dat$sample, keep = dat$keep)) keep
sample FALSE TRUE
3684F1 598 2439
3684F2 1081 7184
3684F3 425 4714
3684F4 795 5266
3684F5 593 5643
3684F6 615 4637
3773F1 7053 8838
3773F10 893 4581
3773F11 1526 4277
3773F12 1495 6245
3773F2 2726 2797
3773F3 1356 3941
3773F4 1585 4410
3773F5 1437 4527
3773F6 1314 5291
3773F7 1668 5127
3773F8 752 2217
3773F9 1031 5297ggplot(as.data.frame(tbl), aes(x = sample, y = Freq, fill = keep)) +
geom_col(position = "stack") + theme_bw() +
theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
labs(x = "", y = "Number of cells")
ggplot(as.data.frame(tbl), aes(x = sample, y = Freq, fill = keep)) +
geom_col(position = "fill") + theme_bw() +
theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
labs(x = "", y = "Fraction of cells")
dat <- dat[, which(dat$keep)]
dim(dat)[1] 20231 87431tbl <- table(dat[[sid]], dat[[stid]])
draw(Heatmap(tbl, col = circlize::colorRamp2(c(0, max(tbl)), c("white", "steelblue")),
cluster_columns = FALSE, cluster_rows = FALSE,
row_names_gp = gpar(fontsize = 8),
name = "Number of\ncells per\nsample",
rect_gp = gpar(col = "grey", lwd = 1),
cell_fun = function(j, i, x, y, width, height, fill) {
grid.text(sprintf("%d", tbl[i, j]), x, y, gp = gpar(fontsize = 8))
}))
Recalculate QC metrics for genes based on remaining cells
dat <- scuttle::addPerFeatureQC(dat, assay.type = "counts")
rowData(dat)$intronic_fraction <-
rowSums(assay(dat, "unspliced")) / rowSums(assay(dat, "total"))Most highly expressed genes
scater::plotHighestExprs(dat, n = 25, colour_cells_by = stid,
exprs_values = "counts")
Calculate size factors and logcounts
size_factor_threshold <- 0.1
set.seed(123)
clusters <- scran::quickCluster(
dat, method = "igraph", min.mean = 0.1,
use.ranks = FALSE, BSPARAM = BiocSingular::ExactParam()
)
dat <- scran::computeSumFactors(dat, min.mean = 0.1, cluster = clusters)
summary(sizeFactors(dat)) Min. 1st Qu. Median Mean 3rd Qu. Max.
0.04152 0.54822 0.91279 1.00000 1.32102 6.98682 ggplot(data.frame(libSize = dat$sum,
sizeFactor = sizeFactors(dat),
group = dat$sample,
stringsAsFactors = FALSE),
aes(x = libSize, y = sizeFactor, color = group)) +
geom_point(alpha = 0.5) + theme_bw() +
geom_hline(yintercept = size_factor_threshold, linetype = "dashed") +
xlab("Total UMI count") + ylab("Size factor")
dat <- dat[, sizeFactors(dat) > size_factor_threshold]
dim(dat)[1] 20231 85628dat <- logNormCounts(dat)Cell cycle assignment
## tricycle
set.seed(123)
dat <- project_cycle_space(dat, species = "human", gname.type = "SYMBOL")
dat <- estimate_cycle_position(dat)
dat <- estimate_Schwabe_stage(dat, gname.type = "SYMBOL", species = "human",
tolerance = 0)
plotReducedDim(dat, "tricycleEmbedding", point_size = 0.5,
colour_by = "CCStage")
tbl <- table(dat$timepoint, dat$CCStage)
draw(Heatmap(tbl, col = circlize::colorRamp2(c(0, max(tbl)), c("white", "steelblue")),
cluster_columns = FALSE, cluster_rows = FALSE,
row_names_gp = gpar(fontsize = 8),
name = "Number of\ncells",
rect_gp = gpar(col = "grey", lwd = 1),
cell_fun = function(j, i, x, y, width, height, fill) {
grid.text(sprintf("%d", tbl[i, j]), x, y, gp = gpar(fontsize = 8))
}))
## Seurat (using genes from HECA)
s_genes <- read.delim(s_genes)[, 1]
g2m_genes <- read.delim(g2m_genes)[, 1]
seurat <- CreateSeuratObject(counts = assay(dat, "counts"))
seurat <- NormalizeData(seurat, verbose = FALSE)
seurat <- FindVariableFeatures(seurat, selection.method = "vst", verbose = FALSE)
seurat <- ScaleData(seurat, features = rownames(seurat), verbose = FALSE)
seurat <- CellCycleScoring(seurat, s.features = s_genes, g2m.features = g2m_genes,
set.ident = TRUE)
stopifnot(all(rownames(seurat@meta.data) == colnames(dat)))
dat$SeuratCC <- seurat@meta.data$Phase
tbl <- table(dat$timepoint, dat$SeuratCC)
draw(Heatmap(tbl, col = circlize::colorRamp2(c(0, max(tbl)), c("white", "steelblue")),
cluster_columns = FALSE, cluster_rows = FALSE,
row_names_gp = gpar(fontsize = 8),
name = "Number of\ncells",
rect_gp = gpar(col = "grey", lwd = 1),
cell_fun = function(j, i, x, y, width, height, fill) {
grid.text(sprintf("%d", tbl[i, j]), x, y, gp = gpar(fontsize = 8))
}))
Find and mark doublets
set.seed(123)
dat <- scDblFinder(dat, samples = "name", verbose = FALSE,
BPPARAM = MulticoreParam(workers = 8L))
table(dat$scDblFinder.class, dat$name)
B050_post_psg_d10 B050_post_psg_d11 B050_post_psg_d12
singlet 4706 4837 3947
doublet 412 451 324
B050_post_psg_d9_12h B050_post_psg_d9_24h B050_pre_psg_d8
singlet 3662 4846 6539
doublet 277 440 760
B066_post_psg_d10 B066_post_psg_d11 B066_post_psg_d12
singlet 4749 5111 4289
doublet 429 519 348
B066_post_psg_d9 B066_post_psg_d9_12h B066_pre_psg_d8
singlet 4332 6349 2297
doublet 375 738 142
B080_post_psg_d10 B080_post_psg_d11 B080_post_psg_d12
singlet 2099 4227 5680
doublet 118 349 557
B080_post_psg_d9_12h B080_post_psg_d9_24h B080_pre_psg_d8
singlet 4060 4163 2633
doublet 349 358 156Sketching
set.seed(123)
geneVar <- scran::modelGeneVar(dat)
subhvgs <- scran::getTopHVGs(geneVar, n = 2000)
subd <- scater::runPCA(dat, exprs_values = "logcounts", ncomponents = 30,
name = "PCAsub",
subset_row = subhvgs, scale = FALSE)
idxsub <- geosketch(reducedDim(subd, "PCAsub"),
N = min(ncol(subd), max(0.1 * ncol(subd), 150)),
seed = 123)
sketchcells <- colnames(subd)[idxsub]
dat$sketched <- ifelse(colnames(dat) %in% sketchcells, TRUE, FALSE)
table(sketched = dat$sketched)sketched
FALSE TRUE
77066 8562 Find an alternative reduced dimension representation, after sketching
set.seed(123)
## First find a new set of HVGs
dattmp <- dat[, sketchcells]
stopifnot(all(rownames(dattmp) == rownames(dat)))
rowData(dat)$scran_geneVar_geosketch <- scran::modelGeneVar(dattmp)
hvgs <- scran::getTopHVGs(rowData(dat)$scran_geneVar_geosketch, n = 3000)
metadata(dat)$hvg.scran.genevar.geosketch <- hvgs
head(
as.data.frame(
rowData(dat)$scran_geneVar_geosketch[
order(rowData(dat)$scran_geneVar_geosketch$bio,
decreasing = TRUE),
c("mean", "total", "tech", "bio", "p.value", "FDR")]
),
25
) mean total tech bio p.value FDR
IFI27 4.424344 7.167424 1.492293 5.675131 2.001221e-64 3.574380e-60
MMP7 4.581938 6.970558 1.492493 5.478065 3.697836e-60 3.302352e-56
CCL20 2.261824 6.167145 1.470837 4.696308 4.922014e-46 2.930403e-42
ISG15 2.580489 5.891723 1.510108 4.381614 2.269297e-38 1.013298e-34
IFI6 2.880933 5.903557 1.532789 4.370768 3.964083e-37 1.416050e-33
IGFBP7 2.213738 5.596281 1.463240 4.133041 1.812628e-36 5.395892e-33
IFITM1 2.716437 5.632226 1.522108 4.110119 1.678515e-33 4.282852e-30
CXCL5 3.173366 5.530787 1.540538 3.990249 5.533206e-31 1.235357e-27
CAPS 3.026924 5.330290 1.537656 3.792634 2.803111e-28 4.551487e-25
CXCL1 2.745286 5.212095 1.524355 3.687740 2.802465e-27 3.850372e-24
TFF3 3.065632 5.201724 1.538530 3.663194 1.746567e-26 1.949715e-23
LCN2 3.534783 5.100962 1.527170 3.573792 1.193569e-25 1.254020e-22
SAA2 1.850499 4.946963 1.391459 3.555504 3.312646e-30 5.916716e-27
GDF15 2.770725 4.706322 1.526227 3.180095 9.916005e-21 7.379574e-18
CXCL8 3.016898 4.556193 1.537494 3.018700 1.288311e-18 8.850203e-16
SAA1 1.716913 4.324135 1.354940 2.969195 9.951645e-23 8.079379e-20
RSPH1 2.175335 4.208529 1.456934 2.751595 2.304293e-17 1.469892e-14
C20orf85 1.306453 3.805729 1.185746 2.619982 4.477304e-23 3.808054e-20
S100A4 2.857985 3.940606 1.531681 2.408925 1.355683e-12 5.044554e-10
VIM 2.388246 3.866303 1.488350 2.377953 6.105068e-13 2.478241e-10
RPS10-NUDT3 2.810523 3.890953 1.528968 2.361985 3.258842e-12 1.164123e-09
CXCL3 1.746871 3.708716 1.363797 2.344919 1.052977e-14 5.224231e-12
SCGB2A1 2.870294 3.865640 1.532276 2.333364 6.439930e-12 2.170257e-09
KRT17 1.321341 3.510942 1.193576 2.317366 3.027459e-18 2.002720e-15
C15orf48 1.087644 3.306294 1.056237 2.250057 1.392162e-21 1.081104e-18## Next run PCA on the sketched cells
set.seed(123)
dattmp <- dat[, sketchcells]
pca <- calculatePCA(dattmp, ncomponents = 30,
subset_row = metadata(dat)$hvg.scran.genevar.geosketch)
projection <- pca
rotation <- attr(pca, "rotation")
centers <- rowMeans(assay(dattmp[metadata(dat)$hvg.scran.genevar.geosketch, ],
"logcounts"))
## Sanity check
stopifnot(rownames(rotation) == names(centers))
tmp <- t(sweep(assay(dattmp[names(centers), ], "logcounts"), MARGIN = 1,
STATS = centers, FUN = "-")) %*% rotation
max(abs(tmp - projection))[1] 1.998401e-13## Project all cells
pca <- t(sweep(assay(dat[names(centers), ], "logcounts"), MARGIN = 1,
STATS = centers, FUN = "-")) %*% rotationWarning in asMethod(object): sparse->dense coercion: allocating vector of size
1.9 GiBstopifnot(all(rownames(pca) == colnames(dat)))
reducedDim(dat, "PCA_sketch") <- as.matrix(pca)Dimension reduction
Apply tSNE and UMAP to sketched cells and project the others.
## TSNE
set.seed(123)
tsne <- snifter::fitsne(
reducedDim(dat[, sketchcells], "PCA_sketch"),
n_jobs = 8L, random_state = 123L,
perplexity = 30, n_components = 2L,
pca = FALSE)
proj <- snifter::project(
x = tsne,
new = reducedDim(dat, "PCA_sketch"),
old = reducedDim(dat[, sketchcells], "PCA_sketch"))
rownames(proj) <- colnames(dat)
reducedDim(dat, "TSNE_sketch") <- proj
## UMAP
set.seed(123)
umap <- uwot::umap(
X = reducedDim(dat[, sketchcells], "PCA_sketch"),
n_components = 2, n_neighbors = 50,
pca = NULL, scale = FALSE,
ret_model = TRUE, n_threads = 8L)
proj <- uwot::umap_transform(
X = reducedDim(dat, "PCA_sketch"),
model = umap,
n_threads = 8L)
rownames(proj) <- colnames(dat)
reducedDim(dat, "UMAP_sketch") <- projShuffle cells (to avoid plotting in layers by sample later)
set.seed(123)
dat <- dat[, sample(seq_len(ncol(dat)), ncol(dat), replace = FALSE)]Clustering
## Helper function
## method should be "louvain" or "leiden"
clusterfcn <- function(sce, dimred = "fastMNNcorr", method = "leiden",
resolutions = c(0.05, 0.1, 0.2, 0.4, 0.6, 1.0),
sketched_cells = NULL) {
if (is.null(sketched_cells)) {
for (r in resolutions) {
set.seed(123)
mat <- reducedDim(sce, dimred)
clst <- bluster::clusterRows(
mat,
BLUSPARAM = bluster::SNNGraphParam(
k = 10, type = "rank",
num.threads = 1L,
cluster.fun = method,
cluster.args = list(resolution = r)))
colData(sce)[[paste0(method, "_res", r)]] <- as.character(clst)
}
} else {
train <- reducedDim(sce[, sketched_cells], dimred)
test <- reducedDim(sce[, setdiff(colnames(sce), sketched_cells)], dimred)
for (r in resolutions) {
set.seed(123)
df <- data.frame(cell = colnames(sce), cluster = NA_character_)
clst <- as.character(bluster::clusterRows(
train,
BLUSPARAM = bluster::SNNGraphParam(
k = 10, type = "rank",
num.threads = 1L,
cluster.fun = method,
cluster.args = list(resolution = r))))
preds <- class::knn(train = train, test = test,
cl = clst, k = 1) |> as.character()
df$cluster[match(rownames(train), df$cell)] <- clst
df$cluster[match(rownames(test), df$cell)] <- preds
stopifnot(colnames(sce) == df$cell)
colData(sce)[[paste0(method, "_res", r)]] <- df$cluster
}
}
sce
}
dat <- clusterfcn(dat, dimred = "PCA_sketch", method = "leiden",
resolutions = c(0.125),
sketched_cells = sketchcells)Annotate clusters
clust_annot <- stack(
list("pre-ciliated" = 7,
"ciliated stress responsive" = 24,
"ciliated wound responsive 12h" = c(23, 15, 27),
"ciliated wound responsive 24h" = 6,
"ciliated immunomodulatory" = c(25, 8),
"ciliated cycling" = 26,
"ciliated" = c(1, 5, 3, 12, 16),
"cycling" = 2,
"stress responsive" = 17,
"wound responsive 12h" = c(19, 20),
"immunomodulatory" = 9,
"wound responsive 24h" = c(10, 11, 14, 22),
"OxPhos high" = 21,
"secretory late" = 18,
"secretory early" = 4,
"hormone responsive" = 13)
)
dat$leiden_res0.125_annot <-
clust_annot$ind[match(dat$leiden_res0.125,
as.character(clust_annot$values))]
table(dat$leiden_res0.125_annot)
pre-ciliated ciliated stress responsive
1046 500
ciliated wound responsive 12h ciliated wound responsive 24h
6598 3030
ciliated immunomodulatory ciliated cycling
5299 821
ciliated cycling
9757 19009
stress responsive wound responsive 12h
398 8645
immunomodulatory wound responsive 24h
3463 16592
OxPhos high secretory late
836 1634
secretory early hormone responsive
2289 5711 leiden_res0.125_annot_cols <- c(
"secretory late" = "firebrick",
"secretory early" = "firebrick1",
"hormone responsive" = "orchid3",
"cycling" = "forestgreen",
"stress responsive" = "grey20",
"wound responsive 12h" = "navajowhite",
"wound responsive 24h" = "goldenrod1",
"OxPhos high" = "coral4",
"immunomodulatory" = "#FF6347",
"pre-ciliated" = "lightskyblue",
"ciliated" = "royalblue",
"ciliated stress responsive" = "navy",
"ciliated wound responsive 12h" = "turquoise",
"ciliated wound responsive 24h" = "turquoise4",
"ciliated immunomodulatory" = "tan1",
"ciliated cycling" = "olivedrab2"
)
atlas_leiden_res0.4_annot_cols <- c(
"pre-ciliated" = "lightskyblue",
"ciliated" = "royalblue",
"secretory early" = "#C8AEE2",
"secretory mid" = "#7D4DBA",
"secretory late" = "#3A1260",
"luminal early" = "gold1",
"luminal late" = "darkorange",
"transcriptionally active" = "turquoise4",
"MUC5B+" = "salmon",
"KRT5+" = "firebrick",
"remodeling" = "#4CAF50",
"cycling" = "mediumaquamarine",
"hormone responsive early" = "hotpink2",
"hormone responsive late" = "hotpink4"
) Dot plots of clusters and marker genes of interest
col = circlize::colorRamp2(
breaks = seq(0, 7.5, length.out = 64),
colors = colorRampPalette(c("grey90", "purple"))(64))
markers <- intersect(
c("MKI67", "MMP7", "STMN1", "PKM", "ALDH1A1", "VIM", "LGR5", "LEFTY1",
"IL6", "ENPP3", "WNT7A", "IL32", "SAA1", "SAA2", "CLU", "MUC5B", "LCN2",
"LY6E", "S100A6", "PLAUR", "PLAU", "SFN", "TNFRSF12A", "IL23A", "KRT17",
"KRT5", "ANXA1", "CXCL8", "SCGB2A1", "SBGB2A2", "SCGB1D2", "IL6ST",
"PAEP", "FGF7", "PGR", "ESR1", "CCNO", "FOXJ1", "PIFO", "RSPH1", "DKK1"),
rownames(dat))
plotDots(dat, features = markers,
group = "leiden_res0.125_annot",
swap_rownames = "gene_name", zlim = c(0, 7.5),
colour = col(seq(0, 7.5, length.out = 64))) +
scale_size(limits=c(0.2, 1)) +
labs(x = "") +
theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5))
Marker genes from HECA paper (https://www.nature.com/articles/s41588-024-01873-w, Fig 2a).
markers <- intersect(
c("PAX8", "EPCAM", "FOXA2", "MMP7", "SOX9", "AXIN2", "TRH", "ALDH1A1",
"CDH2", "KLK11", "SLC7A11", "DKK1", "PHLDA1", "KMO", "IHH", "EMID1",
"HPRT1", "SUFU", "OPRK1", "CBR3", "S100P", "SCGB2A2", "ABCG1", "PAEP",
"DPP4", "GPX3", "PTPRR", "FGF7", "FXYD2", "IL32", "WNT7A", "LGR5",
"MMP7", "TNF", "VTCN1", "CLDN22", "SULT1E1", "LEFTY1", "PTGS1", "IL6",
"CCNO", "CDC20B", "MUC12", "PIFO", "FOXJ1", "TP73", "BPIFB1", "TFF3",
"MUC5B", "SAA1", "TP63", "SNCG", "KRT5", "PGR", "AR", "ESR1"),
rownames(dat))
plotDots(dat, features = markers,
group = "leiden_res0.125_annot",
swap_rownames = "gene_name", zlim = c(0, 7.5),
colour = col(seq(0, 7.5, length.out = 64))) +
scale_size(limits=c(0.2, 1)) +
labs(x = "") +
theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5))
Plot dimension reductions
timepoint_cols <- c(
"Horm Withdr 48h" = "grey40",
"Post-breakdown 12h" = "#f3a9c0",
"Post-breakdown 24h" = "#e27897",
"Post-breakdown 48h" = "#cb3b63",
"E2-diff 24h" = "#debaf5",
"E2-diff 48h" = "#BE9FD2"
)
donor_cols <- c(
"B050" = "#FF9966",
"B066" = "slateblue",
"B080" = "darkgreen"
)
seurat_cols <- c(
"S" = "gold",
"G1" = "grey80",
"G2M" = "brown"
)
for (vrb in c("scDblFinder.class", "leiden_res0.125")) {
print(plotReducedDim(dat, dimred = "UMAP_sketch",
colour_by = vrb, point_size = 0.5) +
guides(color = guide_legend(override.aes = list(size = 2, alpha = 1))) +
theme_void())
}
print(plotReducedDim(dat, dimred = "UMAP_sketch",
colour_by = "leiden_res0.125_annot", point_size = 0.5) +
scale_color_manual(values = leiden_res0.125_annot_cols, name = "") +
guides(color = guide_legend(override.aes = list(size = 2, alpha = 1))) +
theme_void())Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.
print(plotReducedDim(dat, dimred = "UMAP_sketch",
colour_by = "timepoint", point_size = 0.5) +
scale_color_manual(values = timepoint_cols, name = "") +
guides(color = guide_legend(override.aes = list(size = 2, alpha = 1))) +
theme_void())Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.
print(plotReducedDim(dat, dimred = "UMAP_sketch",
colour_by = "donorCol", point_size = 0.5) +
scale_color_manual(values = donor_cols, name = "") +
guides(color = guide_legend(override.aes = list(size = 2, alpha = 1))) +
theme_void())Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.
print(plotReducedDim(dat, dimred = "UMAP_sketch",
colour_by = "SeuratCC", point_size = 0.5) +
scale_color_manual(values = seurat_cols, name = "") +
guides(color = guide_legend(override.aes = list(size = 2, alpha = 1))) +
theme_void())Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.
for (vrb in c("subsets_MT_percent", "detected", "intronic_fraction",
"CCNO", "FOXJ1", "MKI67")) {
print(plotReducedDim(dat, dimred = "UMAP_sketch",
colour_by = vrb, point_size = 0.5) +
theme_void())
}
Transfer labels from tissue atlas
set.seed(123)
atlas <- readRDS(tissue_atlas_sce)
ref_summed <- aggregateAcrossCells(atlas, id = atlas$leiden_res0.4_annot,
use.assay.type = "logcounts",
statistics = "mean")
pred <- SingleR(test = dat, ref = ref_summed,
labels = ref_summed$leiden_res0.4_annot)
plotScoreHeatmap(pred)
stopifnot(all(rownames(pred) == colnames(dat)))
dat$pred_from_invivo <- pred$labels
dat$pred_from_invivo_pruned <- pred$pruned.labels
dat$pred_from_invivo_delta <- pred$delta.nextVisualize predicted labels
print(plotReducedDim(dat, dimred = "UMAP_sketch",
colour_by = "pred_from_invivo", point_size = 0.5) +
scale_color_manual(values = atlas_leiden_res0.4_annot_cols, name = "") +
guides(color = guide_legend(override.aes = list(size = 2, alpha = 1))))Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.
Plot cluster abundance over time
for (col in c("leiden_res0.125_annot", "pred_from_invivo")) {
df <- data.frame(colData(dat)[, c("timepoint", col)]) |>
# mutate(timepoint = factor(timepoint,
# levels = c("d8", "d9_12h", "d9_24h",
# "d10", "d11", "d12"))) |>
group_by(timepoint, .data[[col]]) |>
tally() |>
group_by(timepoint) |>
mutate(perc = n / sum(n) * 100) |>
group_by(.data[[col]]) |>
mutate(scaledperc = perc / sum(perc)) |>
ungroup()
# by donor
dfdonor <- data.frame(colData(dat)[, c("timepoint", "donor", col)]) |>
group_by(timepoint, donor, .data[[col]]) |>
tally() |>
group_by(timepoint, donor) |>
mutate(perc = n / sum(n) * 100) |>
group_by(.data[[col]]) |>
mutate(scaledperc = perc / sum(perc)) |>
ungroup()
if (col == "leiden_res0.125_annot") {
df$cluster_group <- ifelse(grepl("ciliated", df[[col]]),
"ciliated", "non-ciliated")
}
gg <- ggplot(df, aes(x = timepoint, y = perc, color = .data[[col]])) +
geom_point() +
geom_line(aes(x = as.numeric(timepoint)), linewidth = 1) +
theme_cowplot() +
labs(x = "Timepoint", y = "Relative abundance (%)",
title = col) +
theme(legend.position = "bottom")
if (col == "leiden_res0.125_annot") {
gg <- gg + facet_wrap(~ cluster_group, ncol = 1) +
scale_color_manual(values = leiden_res0.125_annot_cols, name = "")
} else {
gg <- gg +
scale_color_manual(values = atlas_leiden_res0.4_annot_cols, name = "")
}
print(gg)
## One panel per cluster
gg <- ggplot(df, aes(x = timepoint, y = perc, color = .data[[col]])) +
geom_point() +
geom_line(aes(x = as.numeric(timepoint)), linewidth = 1) +
theme_cowplot() +
labs(x = "Timepoint", y = "Relative abundance (%)",
title = col) +
facet_wrap(~ .data[[col]], scales = "free_y") +
theme(legend.position = "none",
axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5))
if (col == "leiden_res0.125_annot") {
gg <- gg + scale_color_manual(values = leiden_res0.125_annot_cols, name = "")
} else {
gg <- gg +
scale_color_manual(values = atlas_leiden_res0.4_annot_cols, name = "")
}
print(gg)
## One panel per cluster, add lines for donors
gg <- ggplot(bind_rows(df |> mutate(donor = "overall"), dfdonor),
aes(x = timepoint, y = perc, group = donor, alpha = donor == "overall",
linewidth = donor, color = .data[[col]])) +
geom_point() +
geom_line(aes(x = as.numeric(timepoint)), linewidth = 1) +
scale_alpha_manual(values = c(`TRUE` = 1, `FALSE` = 0.25)) +
scale_linewidth_manual(values = c(`TRUE` = 2, `FALSE` = 0.75)) +
theme_cowplot() +
labs(x = "Timepoint", y = "Relative abundance (%)",
title = col) +
facet_wrap(~ .data[[col]], scales = "free_y") +
theme(legend.position = "none",
axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5))
if (col == "leiden_res0.125_annot") {
gg <- gg + scale_color_manual(values = leiden_res0.125_annot_cols, name = "")
} else {
gg <- gg +
scale_color_manual(values = atlas_leiden_res0.4_annot_cols, name = "")
}
print(gg)
## Bar plot
gg <- ggplot(df, aes(x = timepoint, y = perc, fill = .data[[col]])) +
geom_col() +
theme_cowplot() +
labs(x = "Timepoint", y = "Relative abundance (%)", title = col) +
coord_flip() +
theme(legend.position = "bottom")
if (col == "leiden_res0.125_annot") {
gg <- gg + scale_fill_manual(values = leiden_res0.125_annot_cols, name = "")
} else {
gg <- gg +
scale_fill_manual(values = atlas_leiden_res0.4_annot_cols, name = "")
}
print(gg)
}
Find marker genes
First, we compare all clusters in a pairwise manner. The comparison is only performed for genes with a non-zero count in at least 0.05% of the cells (corresponding to 43 of the 85628 cells in the data set). These gene IDs are also written to a file, for potential use as the ‘universe’ in over-representation analyses.
frac_detected <- rowMeans(assay(dat, "counts") > 0)
genes_to_test <- rownames(dat)[frac_detected >= (0.05 / 100)]
write.table(genes_to_test, file = "organoid_timecourse_gene_universe.txt",
row.names = FALSE, col.names = FALSE, quote = FALSE, sep = "\t")
clc <- "leiden_res0.125_annot"
ptts <- pairwiseTTests(dat, groups = dat[[clc]],
direction = "up",
subset.row = genes_to_test)
# markers wrt all other clusters
markers <- combineMarkers(ptts$statistics, ptts$pairs, pval.type = "all")
# markers within lineage
keep_ciliated <- which(grepl("ciliated", ptts$pairs$first) &
grepl("ciliated", ptts$pairs$second))
markers_ciliated <- combineMarkers(ptts$statistics[keep_ciliated],
ptts$pairs[keep_ciliated, ],
pval.type = "all")
keep_nonciliated <- which(!grepl("ciliated", ptts$pairs$first) &
!grepl("ciliated", ptts$pairs$second))
markers_nonciliated <- combineMarkers(ptts$statistics[keep_nonciliated],
ptts$pairs[keep_nonciliated, ],
pval.type = "all")
## Write results to files
ttestdir <- paste0("pairwise_ttests_", clc)
if (!file.exists(ttestdir)) {
dir.create(ttestdir, recursive = TRUE)
}
for (i in seq_len(nrow(ptts$pairs))) {
df <- cbind(gene = rownames(ptts$statistics[[i]]),
ptts$statistics[[i]])
df <- df[order(df$p.value, -df$logFC), ]
filename <- paste0(clc, "_cluster_", ptts$pairs$first[i],
"_vs_", ptts$pairs$second[i], ".txt")
write.table(df, file = file.path(ttestdir, filename), quote = FALSE,
sep = "\t", row.names = FALSE, col.names = TRUE)
}
markerdir <- paste0("marker_genes_all_", clc)
if (!file.exists(markerdir)) {
dir.create(markerdir, recursive = TRUE)
}
for (i in seq_along(markers)) {
df <- cbind(gene = rownames(markers[[i]]), markers[[i]])
filename <- paste0(clc, "_cluster_", names(markers)[i], "_markers.txt")
write.table(df, file = file.path(markerdir, filename), quote = FALSE,
sep = "\t", row.names = FALSE, col.names = TRUE)
}
markerdirlin <- paste0("marker_genes_within_lineage_all_", clc)
if (!file.exists(markerdirlin)) {
dir.create(markerdirlin, recursive = TRUE)
}
for (i in seq_along(markers_ciliated)) {
df <- cbind(gene = rownames(markers_ciliated[[i]]),
markers_ciliated[[i]])
filename <- paste0(clc, "_cluster_", names(markers_ciliated)[i],
"_markers_within_lineage.txt")
write.table(df, file = file.path(markerdirlin, filename), quote = FALSE,
sep = "\t", row.names = FALSE, col.names = TRUE)
}
for (i in seq_along(markers_nonciliated)) {
df <- cbind(gene = rownames(markers_nonciliated[[i]]),
markers_nonciliated[[i]])
filename <- paste0(clc, "_cluster_", names(markers_nonciliated)[i],
"_markers_within_lineage.txt")
write.table(df, file = file.path(markerdirlin, filename), quote = FALSE,
sep = "\t", row.names = FALSE, col.names = TRUE)
}Save object
saveRDS(dat, file = "organoid_scrnaseq_timecourse.rds")Session info
sessionInfo()R version 4.4.2 (2024-10-31)
Platform: x86_64-pc-linux-gnu
Running under: Rocky Linux 8.10 (Green Obsidian)
Matrix products: default
BLAS/LAPACK: FlexiBLAS OPENBLAS; LAPACK version 3.12.0
locale:
[1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
[3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
[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: Europe/Zurich
tzcode source: system (glibc)
attached base packages:
[1] parallel grid stats4 stats graphics grDevices utils
[8] datasets methods base
other attached packages:
[1] SingleR_2.8.0 swissknife_0.43
[3] DT_0.33 igraph_2.1.4
[5] BiocParallel_1.40.2 bluster_1.16.0
[7] org.Hs.eg.db_3.20.0 AnnotationDbi_1.68.0
[9] tricycle_1.14.0 sketchR_1.2.0
[11] Seurat_5.2.1 SeuratObject_5.0.2
[13] sp_2.2-0 scDblFinder_1.20.2
[15] circlize_0.4.16 ComplexHeatmap_2.22.0
[17] cowplot_1.1.3 msigdbr_10.0.2
[19] tidyr_1.3.1 dplyr_1.1.4
[21] scran_1.34.0 scater_1.34.1
[23] ggplot2_3.5.2 scuttle_1.16.0
[25] SingleCellExperiment_1.28.1 SummarizedExperiment_1.36.0
[27] Biobase_2.66.0 GenomicRanges_1.58.0
[29] GenomeInfoDb_1.42.3 IRanges_2.40.1
[31] S4Vectors_0.44.0 BiocGenerics_0.52.0
[33] MatrixGenerics_1.18.1 matrixStats_1.5.0
[35] fishpond_2.12.0
loaded via a namespace (and not attached):
[1] fs_1.6.6 spatstat.sparse_3.1-0
[3] bitops_1.0-9 httr_1.4.7
[5] RColorBrewer_1.1-3 doParallel_1.0.17
[7] tools_4.4.2 sctransform_0.4.1
[9] R6_2.6.1 lazyeval_0.2.2
[11] uwot_0.2.3 GetoptLong_1.0.5
[13] withr_3.0.2 gridExtra_2.3
[15] progressr_0.15.1 cli_3.6.4
[17] Cairo_1.6-2 spatstat.explore_3.4-2
[19] fastDummies_1.7.5 labeling_0.4.3
[21] sass_0.4.10 mvtnorm_1.3-3
[23] spatstat.data_3.1-6 readr_2.1.5
[25] ggridges_0.5.6 pbapply_1.7-2
[27] Rsamtools_2.22.0 dichromat_2.0-0.1
[29] parallelly_1.43.0 limma_3.62.2
[31] RSQLite_2.3.9 circular_0.5-1
[33] generics_0.1.3 shape_1.4.6.1
[35] BiocIO_1.16.0 crosstalk_1.2.1
[37] gtools_3.9.5 ica_1.0-3
[39] spatstat.random_3.3-3 Matrix_1.7-1
[41] ggbeeswarm_0.7.2 abind_1.4-8
[43] lifecycle_1.0.4 yaml_2.3.10
[45] edgeR_4.4.2 SparseArray_1.6.2
[47] Rtsne_0.17 blob_1.2.4
[49] promises_1.3.2 dqrng_0.4.1
[51] crayon_1.5.3 dir.expiry_1.14.0
[53] miniUI_0.1.2 lattice_0.22-6
[55] beachmat_2.22.0 KEGGREST_1.46.0
[57] magick_2.8.6 pillar_1.10.2
[59] knitr_1.50 metapod_1.14.0
[61] boot_1.3-31 rjson_0.2.23
[63] xgboost_1.7.9.1 future.apply_1.11.3
[65] codetools_0.2-20 glue_1.8.0
[67] spatstat.univar_3.1-2 data.table_1.17.0
[69] vctrs_0.6.5 png_0.1-8
[71] spam_2.11-1 gtable_0.3.6
[73] assertthat_0.2.1 cachem_1.1.0
[75] xfun_0.52 S4Arrays_1.6.0
[77] mime_0.13 survival_3.7-0
[79] pheatmap_1.0.12 iterators_1.0.14
[81] statmod_1.5.0 fitdistrplus_1.2-2
[83] ROCR_1.0-11 nlme_3.1-166
[85] usethis_3.1.0 bit64_4.6.0-1
[87] filelock_1.0.3 RcppAnnoy_0.0.22
[89] bslib_0.9.0 irlba_2.3.5.1
[91] vipor_0.4.7 KernSmooth_2.23-24
[93] DBI_1.2.3 colorspace_2.1-1
[95] tidyselect_1.2.1 bit_4.6.0
[97] compiler_4.4.2 curl_6.2.2
[99] BiocNeighbors_2.0.1 basilisk.utils_1.18.0
[101] DelayedArray_0.32.0 plotly_4.10.4
[103] msigdbdf_24.1.0 rtracklayer_1.66.0
[105] scales_1.3.0.9000 lmtest_0.9-40
[107] stringr_1.5.1 digest_0.6.37
[109] goftest_1.2-3 spatstat.utils_3.1-3
[111] rmarkdown_2.29 basilisk_1.18.0
[113] XVector_0.46.0 htmltools_0.5.8.1
[115] pkgconfig_2.0.3 sparseMatrixStats_1.18.0
[117] fastmap_1.2.0 rlang_1.1.6
[119] GlobalOptions_0.1.2 htmlwidgets_1.6.4
[121] UCSC.utils_1.2.0 DelayedMatrixStats_1.28.1
[123] shiny_1.10.0 jquerylib_0.1.4
[125] farver_2.1.2 zoo_1.8-14
[127] jsonlite_2.0.0 BiocSingular_1.22.0
[129] RCurl_1.98-1.17 magrittr_2.0.3
[131] GenomeInfoDbData_1.2.13 dotCall64_1.2
[133] patchwork_1.3.0 Rcpp_1.0.14
[135] ggnewscale_0.5.1 babelgene_22.9
[137] viridis_0.6.5 reticulate_1.42.0
[139] stringi_1.8.7 zlibbioc_1.52.0
[141] MASS_7.3-61 plyr_1.8.9
[143] listenv_0.9.1 ggrepel_0.9.6
[145] deldir_2.0-4 Biostrings_2.74.1
[147] splines_4.4.2 tensor_1.5
[149] hms_1.1.3 locfit_1.5-9.12
[151] spatstat.geom_3.3-6 RcppHNSW_0.6.0
[153] reshape2_1.4.4 ScaledMatrix_1.14.0
[155] XML_3.99-0.18 evaluate_1.0.3
[157] snifter_1.16.0 tzdb_0.5.0
[159] foreach_1.5.2 httpuv_1.6.16
[161] RANN_2.6.2 purrr_1.0.4
[163] polyclip_1.10-7 future_1.40.0
[165] clue_0.3-66 scattermore_1.2
[167] rsvd_1.0.5 xtable_1.8-4
[169] restfulr_0.0.15 svMisc_1.4.3
[171] RSpectra_0.16-2 later_1.4.2
[173] class_7.3-22 viridisLite_0.4.2
[175] tibble_3.2.1 memoise_2.0.1
[177] beeswarm_0.4.0 GenomicAlignments_1.42.0
[179] cluster_2.1.6 globals_0.17.0