Last updated: 2026-04-22

Checks: 6 1

Knit directory: cf-eti-bal-scrna/

This reproducible R Markdown analysis was created with workflowr (version 1.7.0). The Checks tab describes the reproducibility checks that were applied when the results were created. The Past versions tab lists the development history.

R Markdown file: uncommitted changes

The R Markdown file has unstaged changes. To know which version of the R Markdown file created these results, you’ll want to first commit it to the Git repo. If you’re still working on the analysis, you can ignore this warning. When you’re finished, you can run wflow_publish to commit the R Markdown file and build the HTML.

Environment: empty

Great job! The global environment was empty. Objects defined in the global environment can affect the analysis in your R Markdown file in unknown ways. For reproduciblity it’s best to always run the code in an empty environment.

Seed: set.seed(20260422)

The command set.seed(20260422) was run prior to running the code in the R Markdown file. Setting a seed ensures that any results that rely on randomness, e.g. subsampling or permutations, are reproducible.

Session information: recorded

Great job! Recording the operating system, R version, and package versions is critical for reproducibility.

Cache: none

Nice! There were no cached chunks for this analysis, so you can be confident that you successfully produced the results during this run.

File paths: relative

Great job! Using relative paths to the files within your workflowr project makes it easier to run your code on other machines.

Repository version: 8ac1a7b

Great! You are using Git for version control. Tracking code development and connecting the code version to the results is critical for reproducibility.

The results in this page were generated with repository version 8ac1a7b. See the Past versions tab to see a history of the changes made to the R Markdown and HTML files.

Note that you need to be careful to ensure that all relevant files for the analysis have been committed to Git prior to generating the results (you can use wflow_publish or wflow_git_commit). workflowr only checks the R Markdown file, but you know if there are other scripts or data files that it depends on. Below is the status of the Git repository when the results were generated:


Ignored files:
    Ignored:    .DS_Store
    Ignored:    ._.DS_Store
    Ignored:    code/._bal_RecM.integration_clustering.R
    Ignored:    data/._README.md
    Ignored:    data/AM.markers.csv
    Ignored:    data/DE/
    Ignored:    data/Major.markers.csv
    Ignored:    data/RecM.markers.csv
    Ignored:    data/SCEs
    Ignored:    data/cellbender
    Ignored:    data/cytokine
    Ignored:    data/plots
    Ignored:    data/sample_sheets
    Ignored:    data/vireo
    Ignored:    output/plots/

Untracked files:
    Untracked:  code/scVelo.ipynb
    Untracked:  data/annotation/
    Untracked:  data/public_datasets

Unstaged changes:
    Modified:   analysis/01.demux_BP1.Rmd
    Modified:   analysis/01.demux_BP2.Rmd
    Modified:   analysis/01.demux_Pi1.Rmd
    Modified:   analysis/01.demux_capture1.Rmd
    Modified:   analysis/02.qc_BP1.Rmd
    Modified:   analysis/02.qc_BP2.Rmd
    Modified:   analysis/02.qc_Pi1.Rmd
    Modified:   analysis/02.qc_capture1.Rmd
    Modified:   analysis/04.annotation.Rmd
    Modified:   analysis/05.RecM_analysis.Rmd
    Modified:   analysis/06.MDM-PLA2G7_analysis.Rmd
    Modified:   analysis/07.AM_analysis.Rmd
    Modified:   analysis/08.ProlifM_analysis.Rmd
    Modified:   analysis/09.ETI_analysis.Rmd
    Modified:   analysis/10.cytokine_analysis.Rmd
    Modified:   analysis/11.cytokine_analysis_ETI.Rmd

Note that any generated files, e.g. HTML, png, CSS, etc., are not included in this status report because it is ok for generated content to have uncommitted changes.

These are the previous versions of the repository in which changes were made to the R Markdown (analysis/01.demux_Pi1.Rmd) and HTML (docs/01.demux_Pi1.html) files. If you’ve configured a remote Git repository (see ?wflow_git_remote), click on the hyperlinks in the table below to view the files as they were in that past version.

File	Version	Author	Date	Message
Rmd	29ed4f7	tcwong1994	2026-04-22	Initial import of workflowr site content and analysis files
html	29ed4f7	tcwong1994	2026-04-22	Initial import of workflowr site content and analysis files

knitr::opts_chunk$set(warning = FALSE, message = FALSE)
xaringanExtra::use_panelset()

1 Load packages

suppressPackageStartupMessages({
library(DropletUtils)
library(here)
library(ggplot2)
library(Seurat)
library(cowplot)
library(patchwork)
library(scater)
library(dplyr)
library(vcfR)
library(janitor)
library(stringr)
library(forcats)
})
set.seed(1990)

2 Set names for this pool

# Specify batch name
batch_name <- "Pi1"
num_of_captures <- 2

# Specify capture name
capture_names <- c(paste0(batch_name, "-c",1:num_of_captures))
capture_names <- setNames(capture_names, capture_names)

# Assign sample ID to HTO ID
# Manually list sample names matching HTO id (in alphabetical order)
samples <- c("M1C177", #HTO3
             "M1C177C", #HTO6
             "M1C170C", #HTO10
             "M1C176", #HTO14
             "M1C180", #HTO15
             "M1C180D" #HTO16
             )

3 Read CellBender outputs

Ambient RNA and empty droplets were removed by CellBender (see code/cellbender.sh).

captures <- setNames(
  here("data",
       "cellbender",
       capture_names,
       paste0(capture_names,".cellbender_filtered.h5")),
  capture_names)

sce <- read10xCounts(samples = captures, col.names = TRUE)
stopifnot(!anyDuplicated(colnames(sce)))

sce <- splitAltExps(
  sce,
  rowData(sce)$Type,
  "Gene Expression")

# Tidy up colData
sce$Capture <- factor(sce$Sample)
capture_names <- levels(sce$Capture)
capture_names <- setNames(capture_names, capture_names)
sce$Sample <- NULL

# Number of droplets after CellBender
dim(sce)

[1] 36601 55519

4 HTO demultiplexing

HTODemux from the Seurat package was used to demultiplex droplets. Demultiplexing was performed separately on each capture.

4.1 Prepare HTO data

# change DelayedMatrix to dgCMatrix
assay(sce, "counts") <- as(assay(sce, "counts"), "dgCMatrix")

# Read features.csv used in cellranger multi
features <- read.csv(here("data",
                          "sample_sheets",
                          paste0(batch_name,".features.csv")))

# HTO id are always in alphabetical order
features

            id         name read pattern        sequence     feature_type
1  Human_HTO_3  Human_HTO_3   R2  5P(BC) TTCCGCCTCTCTTTG Antibody Capture
2  Human_HTO_6  Human_HTO_6   R2  5P(BC) GGTTGCCAGATGTCA Antibody Capture
3 Human_HTO_10 Human_HTO_10   R2  5P(BC) ATTGACCCGCGTTAG Antibody Capture
4 Human_HTO_14 Human_HTO_14   R2  5P(BC) CTGTATGTCCGATTG Antibody Capture
5 Human_HTO_15 Human_HTO_15   R2  5P(BC) TAAGATTCAGAGCGA Antibody Capture
6 Human_HTO_16 Human_HTO_16   R2  5P(BC) CTCAGTGCATTCTGG Antibody Capture

# Rename assay
is_hto <- grepl("^Human_HTO", rownames(altExp(sce, "Antibody Capture")))
altExp(sce, "HTO") <- altExp(sce, "Antibody Capture")[is_hto, ]
altExp(sce, "Antibody Capture") <- NULL

4.2 HTO Demultiplexing

# load in umi and hto matrix
umis <- counts(sce)
htos <- counts(altExp(sce))

# Select cell barcodes detected by both RNA and HTO
joint.bcs <- intersect(colnames(umis), colnames(htos))

# Subset RNA and HTO counts by joint cell barcodes
umis <- umis[, joint.bcs]
htos <- as.matrix(htos[, joint.bcs])

# Confirm that the HTO have the correct names
rownames(htos)

[1] "Human_HTO_3"  "Human_HTO_6"  "Human_HTO_10" "Human_HTO_14" "Human_HTO_15"
[6] "Human_HTO_16"

# Perform demultiplexing separately on each capture
obj.list <- list()

lappend <- function (lst, ...){
  lst <- c(lst, list(...))
  return(lst)
}

for (cn in capture_names) {
  message(cn)
  umi = as.matrix(umis[, sce$Capture == cn])
  hto = as.matrix(htos[, sce$Capture == cn])
  
  # Setup Seurat object
  hashtag <- CreateSeuratObject(counts = umi)
  
  # Normalize RNA data with log normalization
  hashtag <- NormalizeData(hashtag)
  
  # Find and scale variable features
  hashtag <- FindVariableFeatures(hashtag, selection.method = "mean.var.plot")
  hashtag <- ScaleData(hashtag, features = VariableFeatures(hashtag))
  
  # Add HTO data as a new assay independent from RNA
  hashtag[["HTO"]] <- CreateAssayObject(counts = hto)
  
  # Normalize HTO data, here we use centered log-ratio (CLR) transformation
  hashtag <- NormalizeData(hashtag, assay = "HTO", normalization.method = "CLR")
  
  # Demultiplex cells based on HTO enrichment
  hashtag <- HTODemux(hashtag, assay = "HTO", positive.quantile = 0.99)
  
  # Rename hash.ID column 
  hashtag$hash.ID <- factor(gsub("-","_",hashtag$hash.ID), 
                            levels = c(features$id, "Doublet", "Negative"))
  obj.list <- lappend(obj.list, hashtag)
}

names(obj.list) <- capture_names
hashtag <- NULL

4.3 Add HTO demultiplexing results to SCE

# add hash.ID column to SCE object
hash.ID <- Reduce(c,
                  lapply(obj.list, function(x) 
                    setNames(x@meta.data$hash.ID, colnames(x))))
sce$HTODemux.ID <- hash.ID

# Store HTODemux results to SCE object
sce$HTODemux_result <- bind_rows(lapply(obj.list, function(x) x@meta.data))

# Add sample ID to SCE object
sampleID.HTO <- setNames(c(samples,"Doublet","Negative"), levels(sce$HTODemux.ID))
sampleID.HTO <- sampleID.HTO[sce$HTODemux.ID]
names(sampleID.HTO) <- colnames(sce)
sce$sampleID.HTO <- factor(sampleID.HTO, levels=c(samples,"Doublet","Negative"))

5 Visualize HTO demultiplexing results

Some code in this section are adapted from internal code by Dr. Peter F. Hickey.

# add colour
sce$colours <- S4Vectors::make_zero_col_DFrame(ncol(sce))

hto_colours <- setNames(
  c(palette.colors(nlevels(sce$HTODemux.ID)-2, "Tableau 10"),"#ced4da","#6c757d"),
  levels(sce$HTODemux.ID))
sce$colours$hto_colours <- hto_colours[sce$HTODemux.ID]

sample_colours <- setNames(
  c(palette.colors(nlevels(sce$sampleID.HTO)-2, "Tableau 10"),"#ced4da","#6c757d"),
  levels(sce$sampleID.HTO))
sce$colours$sample_colours <- sample_colours[sce$sampleID.HTO]

capture_colours <- setNames(
  palette.colors(nlevels(sce$Capture), "Accent"),
  levels(sce$Capture))
sce$colours$capture_colours <- capture_colours[sce$Capture]

5.1 Number of singlets

Pi1-c1

table(obj.list$`Pi1-c1`$HTO_classification.global)


 Doublet Negative  Singlet 
    4928     6187    17748

Pi1-c2

table(obj.list$`Pi1-c2`$HTO_classification.global)


 Doublet Negative  Singlet 
    4149     4772    17735

5.2 HTO expression in ridge plot

Pi1-c1

# Group cells based on the max HTO signal
Idents(obj.list$`Pi1-c1`) <- "hash.ID"
RidgePlot(obj.list$`Pi1-c1`, assay = "HTO", 
          features = rownames(obj.list$`Pi1-c1`[["HTO"]]))

Version	Author	Date
29ed4f7	tcwong1994	2026-04-22

Pi1-c2

# Group cells based on the max HTO signal
Idents(obj.list$`Pi1-c2`) <- "hash.ID"
RidgePlot(obj.list$`Pi1-c2`, assay = "HTO", 
          features = rownames(obj.list$`Pi1-c2`[["HTO"]]))

Version	Author	Date
29ed4f7	tcwong1994	2026-04-22

5.3 HTO expression in heatmap

Pi1-c1

HTOHeatmap(obj.list$`Pi1-c1`, assay = "HTO", ncells = 5000) + 
  ggtitle("Pi1-c1") + 
  theme(plot.title = element_text(hjust=0.5, size=10))

Version	Author	Date
29ed4f7	tcwong1994	2026-04-22

Pi1-c2

HTOHeatmap(obj.list$`Pi1-c2`, assay = "HTO", ncells = 5000) + 
  ggtitle("Pi1-c2") + 
  theme(plot.title = element_text(hjust=0.5, size=12))

Version	Author	Date
29ed4f7	tcwong1994	2026-04-22

5.4 tSNE embeddings for HTOs

Pi1-c1

# remove negative cells from the object
hashtag.subset <- subset(obj.list$`Pi1-c1`, idents = "Negative", invert = TRUE)

# calculate a tSNE embedding of the HTO data
DefaultAssay(hashtag.subset) <- "HTO"
hashtag.subset <- ScaleData(hashtag.subset, features = rownames(hashtag.subset))
hashtag.subset <- RunPCA(hashtag.subset, features = rownames(hashtag.subset), approx = FALSE)
hashtag.subset <- RunTSNE(hashtag.subset, dims = 1:length(rownames(hashtag.subset)), 
                          perplexity = 100, check_duplicates = FALSE)

# tSNE plot
DimPlot(hashtag.subset)

Version	Author	Date
29ed4f7	tcwong1994	2026-04-22

Pi1-c2

# remove negative cells from the object
hashtag.subset <- subset(obj.list$`Pi1-c2`, idents = "Negative", invert = TRUE)

# calculate a tSNE embedding of the HTO data
DefaultAssay(hashtag.subset) <- "HTO"
hashtag.subset <- ScaleData(hashtag.subset, features = rownames(hashtag.subset))
hashtag.subset <- RunPCA(hashtag.subset, features = rownames(hashtag.subset), approx = FALSE)
hashtag.subset <- RunTSNE(hashtag.subset, dims = 1:length(rownames(hashtag.subset)), 
                          perplexity = 100, check_duplicates = FALSE)

# tSNE plot
DimPlot(hashtag.subset)

Version	Author	Date
29ed4f7	tcwong1994	2026-04-22

5.5 Number of demultiplexed droplets

df <- as.data.frame(colData(sce))

Named by HTO

# Number of droplets per HTO classification
p1 <- ggplot(df) + 
  geom_bar(
    aes(
      x = factor(hash.ID, levels(sce$HTODemux.ID)),
      fill=hash.ID), 
    position = position_stack(reverse = TRUE)) + 
  geom_text(stat='count', aes(x = hash.ID, label=..count..), hjust=1, size=2.5) +
  coord_flip() +
  scale_fill_manual(values = hto_colours) +
  ylab("Number of droplets") +
  xlab("HTODemux classification") +
  theme_cowplot(font_size = 8) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

p1 + NoLegend() + p1 + facet_grid(~sce$Capture) + plot_layout(widths=c(1, length(capture_names)))

Version	Author	Date
29ed4f7	tcwong1994	2026-04-22

Named by sample ID

# Number of droplets per sample
p2 <- ggplot(df) + 
  geom_bar(
    aes(
      x = factor(sampleID.HTO, levels(sce$sampleID.HTO)),
      fill=sampleID.HTO), 
    position = position_stack(reverse = TRUE)) + 
  geom_text(stat='count', aes(x = sampleID.HTO, label=..count..), hjust=1, size=2.5) +
  coord_flip() +
  scale_fill_manual(values = sample_colours) +
  ylab("Number of droplets") +
  xlab("Sample ID") +
  theme_cowplot(font_size = 8) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

p2 + NoLegend() + p2 + facet_grid(~sce$Capture) + plot_layout(widths=c(1, length(capture_names)))

Version	Author	Date
29ed4f7	tcwong1994	2026-04-22

5.6 Proportion of singlets per capture

p3 <- ggplot(df) + 
  geom_bar(
    aes(
      x = Capture,
      fill= HTODemux_result.HTO_classification.global), 
    position = position_fill(reverse = TRUE)) + 
  coord_flip() +
  scale_fill_manual(values = setNames(c("#6c757d","#ced4da","#cfe1b9"),
                                      c("Negative","Doublet","Singlet"))) +
  ylab("Proportion of droplets") +
  xlab("HTODemux classification") +
  theme_cowplot(font_size = 8) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

p3 +
  guides(fill=guide_legend(title="HTO_classification.global"))

Version	Author	Date
29ed4f7	tcwong1994	2026-04-22

6 Genetic demultiplexing

Some code in this section are adapted from internal code by Dr. Peter F. Hickey.

SNPs were called from each capture using cellsnp-lite v1.2.3. Based on the SNP information, cells were assigned to genetic donors using vireo v0.5.8.

6.1 Read vireo results

vireo_df <- do.call(
  rbind,
  lapply(capture_names, function(cn) {
    message(cn)
    vireo_df <- read.table(
      here("data", "vireo", cn, "donor_ids.tsv"),
      header = TRUE)
    # Rename 'doublet' and 'unassigned' to match the terms used in HTODemux 
    vireo_df$donor_id <- gsub("doublet", "Doublet", vireo_df$donor_id)
    vireo_df$donor_id <- gsub("unassigned", "Negative", vireo_df$donor_id)
    vireo_df$donor_id <- paste0(cn, "_", vireo_df$donor_id)
    captureNumber <- str_sub(cn, start= -1)
    vireo_df$colname <- paste0(captureNumber, "_", vireo_df$cell)
    # reorder columns to matches SCE.
    j <- match(colnames(sce)[sce$Capture == cn], vireo_df$colname)
    stopifnot(!anyNA(j))
    vireo_df <- vireo_df[j, ]
  }))

6.2 Summarise vireo results

knitr::kable(
  tabyl(vireo_df, donor_id) %>%
    adorn_pct_formatting(1),
  caption = "Assignment of droplets to donors using vireo.")

Assignment of droplets to donors using vireo.
donor_id	n	percent
Pi1-c1_Doublet	7531	13.6%
Pi1-c1_Negative	3887	7.0%
Pi1-c1_donor0	2408	4.3%
Pi1-c1_donor1	8359	15.1%
Pi1-c1_donor2	5443	9.8%
Pi1-c1_donor3	1235	2.2%
Pi1-c2_Doublet	4652	8.4%
Pi1-c2_Negative	4571	8.2%
Pi1-c2_donor0	2376	4.3%
Pi1-c2_donor1	8357	15.1%
Pi1-c2_donor2	5331	9.6%
Pi1-c2_donor3	1369	2.5%

6.3 Identify the best match between donors from the combination of captures

# capture 1 is used as base
cn1 <- capture_names[1]

# read vcf
f1 <- here("data","vireo",cn1,"GT_donors.vireo.vcf.gz")
x1 <- read.vcfR(f1, verbose=FALSE)

# create unique ID for each locus in each capture.
y1 <- paste(
  x1@fix[,"CHROM"],
  x1@fix[,"POS"],
  x1@fix[,"REF"],
  x1@fix[,"ALT"],
  sep = "_")

# match donors of every remaining captures to capture 1
f.best_match_df <- data.frame()

for (cn2 in capture_names[2:length(capture_names)]) {
  
  # read vcf
  f2 <- here("data","vireo",cn2,"GT_donors.vireo.vcf.gz")
  x2 <- read.vcfR(f2, verbose=FALSE)
  
  # create unique ID for each locus in each capture.
  y2 <- paste(
    x2@fix[,"CHROM"],
    x2@fix[,"POS"],
    x2@fix[,"REF"],
    x2@fix[,"ALT"],
    sep = "_")
  
  # only keep the loci in common between the 2 captures.
  i1 <- na.omit(match(y2, y1))
  i2 <- na.omit(match(y1, y2))
  
  # construct genotype matrix at common loci from the 2 captures.
  donor_names <- colnames(x1@gt)[-1]
  
  g1 <- apply(
    x1@gt[i1, donor_names],
    2,
    function(x) sapply(strsplit(x, ":"), `[[`, 1))
  g2 <- apply(
    x2@gt[i2, donor_names],
    2,
    function(x) sapply(strsplit(x, ":"), `[[`, 1))
  
  # count number of genotype matches between pairs of donors (one from each 
  # capture) and convert to a proportion.
  z <- matrix(
    NA_real_,
    nrow = length(donor_names),
    ncol = length(donor_names),
    dimnames = list(donor_names, donor_names))
  
  for (i in rownames(z)) {
    for (j in colnames(z)) {
      z[i, j] <- sum(g1[, i] == g2[, j]) / nrow(g1)
    }
  }
  rownames(z) <- paste(cn1, rownames(z),sep="_")
  colnames(z) <- paste(cn2, colnames(z),sep="_")
  
  # look for the best match between donors
  best_match_df <- data.frame(
    cn1 = rownames(z),
    cn2 = apply(
      z, 1,
      function(x) colnames(z)[which.max(x)]),
    check.names = FALSE)
  
  knitr::kable(
    dplyr::select(best_match_df, everything()),
    caption = paste0("Best match between donors between ", 
                     cn1, " and ", cn2, "."),
    row.names = FALSE)
  
  stopifnot(identical(colnames(sce), vireo_df$colname))
  
  # visualize the best match in a heatmap
  pheatmap::pheatmap(
    z,
    color = viridisLite::inferno(101),
    cluster_rows = FALSE,
    cluster_cols = FALSE,
    main = paste0("Proportion of matching genotypes between ", 
                  cn1, " and ", cn2, "."))
  
  # add matching results to the final best match data frame
  if (nrow(f.best_match_df)==0) {
    f.best_match_df <- best_match_df
  } else {
    f.best_match_df <- f.best_match_df %>% 
      mutate(c=best_match_df$cn2)
    colnames(f.best_match_df) <- paste0("cn",1:ncol(f.best_match_df))
  }
}

Version	Author	Date
29ed4f7	tcwong1994	2026-04-22

# add genetic donor name
f.best_match_df$genetic_donor <- paste0("donor_", LETTERS[1:length(donor_names)])


knitr::kable(
  dplyr::select(f.best_match_df, genetic_donor, everything()),
  caption = "Best match between donors from the two captures.",
  row.names = FALSE)

Best match between donors from the two captures.
genetic_donor	cn1	cn2
donor_A	Pi1-c1_donor0	Pi1-c2_donor0
donor_B	Pi1-c1_donor1	Pi1-c2_donor1
donor_C	Pi1-c1_donor2	Pi1-c2_donor2
donor_D	Pi1-c1_donor3	Pi1-c2_donor3

# transform to a long data frame 
best_match_long_df <- tidyr::pivot_longer(
  data=f.best_match_df,
  cols=starts_with("cn")) %>%
  dplyr::select(genetic_donor, value)

6.4 Add genetic demultiplexing results to SCE

# add genetic donor to SCE object
sce$vireo <- DataFrame(
  vireo_df[, setdiff(colnames(vireo_df), c("cell", "colname"))])

sce$genetic_donor <- left_join(
  vireo_df,
  best_match_long_df, 
  by = c("donor_id" = "value")) %>%
  mutate(
    genetic_donor = factor(
      case_when(
        is.na(genetic_donor) ~ sub(paste0(batch_name,"-c[0-9]_"), "", donor_id),
        TRUE ~ genetic_donor),
      levels = c(paste0("donor_", LETTERS[1:length(donor_names)]), "Doublet", "Negative"))) %>%
  pull(genetic_donor)

# generate a table of HTO to genetic droplets
tb <- as.data.frame.matrix(
  table(as.data.frame(colData(sce)[, c("sampleID.HTO", "genetic_donor")]))) %>%
  dplyr::select(1:length(donor_names)) %>% # discard columns of Doublets and Negatives
  slice_head(n=length(samples)) # discard rows of Doublets and Negatives

# match donor to sample ID according to the best match with HTO
sampleID.genetics <- setNames(c(rownames(tb)[apply(tb, MARGIN = 2, which.max)],
                                "Doublet","Negative"),
                              c(colnames(tb),
                                "Doublet","Negative"))
sampleID.genetics <- sampleID.genetics[sce$genetic_donor]
names(sampleID.genetics) <- colnames(sce)
sce$sampleID.genetics <- factor(sampleID.genetics, 
                                levels=c(samples,"Doublet","Negative"))

7 Visualize genetic demultiplexing results

7.1 Table of number of droplets

janitor::tabyl(
  as.data.frame(colData(sce)[, c("HTODemux.ID", "genetic_donor")]),
  HTODemux.ID,
  genetic_donor) |>
  adorn_title(placement = "combined") |>
  adorn_totals("both") |>
  knitr::kable(
    caption = "Number of droplets assigned to each `HTO`/`Genetic donor` combination.")

Number of droplets assigned to each `HTO`/`Genetic donor` combination.
HTODemux.ID/genetic_donor	donor_A	donor_B	donor_C	donor_D	Doublet	Negative	Total
Human_HTO_3	45	60	1244	14	1093	1429	3885
Human_HTO_6	15	15	7532	7	413	452	8434
Human_HTO_10	3160	31	11	5	690	724	4621
Human_HTO_14	22	25	12	2034	431	776	3300
Human_HTO_15	37	8650	22	6	641	759	10115
Human_HTO_16	17	3867	14	7	399	824	5128
Doublet	533	2310	1362	288	4267	317	9077
Negative	955	1758	577	243	4249	3177	10959
Total	4784	16716	10774	2604	12183	8458	55519

7.2 Number of droplets (named by HTO and genetic donor)

# add colour
genetic_donor_colours <- setNames(
  c(palette.colors(nlevels(sce$genetic_donor)-2, "Set3"), "#d4a276","#9c6644"),
  levels(sce$genetic_donor))
sce$colours$genetic_donor_colours <- genetic_donor_colours[sce$genetic_donor]

# number of droplets assigned by HTO method
p1 <- ggcells(sce) + 
  geom_bar(aes(x = HTODemux.ID, fill = HTODemux.ID)) + 
  geom_text(stat='count', aes(x = HTODemux.ID, label=..count..), 
            hjust=1, size=2) +
  coord_flip() + 
  ylab("Number of droplets") + 
  theme_cowplot(font_size = 8) + 
  scale_y_continuous(breaks=seq(0,16000,4000),limits=c(0,18000)) +
  scale_fill_manual(values = sce$colours$hto_colours) +
  theme(axis.text.x=element_text(angle=45, hjust=1))

p1.facet <- ggcells(sce) + 
  geom_bar(aes(x = HTODemux.ID, fill = HTODemux.ID)) + 
  geom_text(stat='count', aes(x = HTODemux.ID, label=..count..), 
            hjust=1, size=2) +
  #ggtitle("By genetics") +
  ylab("Number of droplets") + 
  facet_grid(~Capture, scales = "fixed", space = "fixed") + 
  theme_cowplot(font_size = 8) + 
  scale_y_continuous(breaks=seq(0,8000,2000), limits=c(0,10000)) +
  scale_fill_manual(values = sce$colours$hto_colours) +
  theme(axis.title.y = element_blank(),
        axis.text.x=element_text(angle=45, hjust=1)) +
  coord_flip()

# number of droplets assigned by genetic method
p2 <- ggcells(sce) + 
  geom_bar(aes(x = genetic_donor, fill = genetic_donor)) + 
  geom_text(stat='count', aes(x = genetic_donor, label=..count..,),
            hjust=1, size=2) +
  coord_flip() + 
  ylab("Number of droplets") + 
  theme_cowplot(font_size = 8) + 
  scale_y_continuous(breaks=seq(0,16000,4000),limits=c(0,18000)) +
  scale_fill_manual(values = sce$colours$genetic_donor_colours) +
  theme(axis.text.x=element_text(angle=45, hjust=1))

p2.facet <- ggcells(sce) + 
  geom_bar(aes(x = genetic_donor, fill = genetic_donor)) + 
  geom_text(stat='count', aes(x = genetic_donor, label=..count..), 
            hjust=1, size=2) +
  #ggtitle("By genetics") +
  ylab("Number of droplets") + 
  facet_grid(~Capture, scales = "fixed", space = "fixed") + 
  theme_cowplot(font_size = 8) + 
  scale_y_continuous(breaks=seq(0,8000,2000), limits=c(0,10000)) +
  scale_fill_manual(values = sce$colours$genetic_donor_colours) +
  theme(axis.title.y = element_blank(),
        axis.text.x=element_text(angle=45, hjust=1)) +
  coord_flip()

(p1 + p1.facet+plot_layout(width=c(1,2))) /
  (p2 + p2.facet+plot_layout(width=c(1,2))) +
  plot_layout(guides="collect")

Version	Author	Date
29ed4f7	tcwong1994	2026-04-22

7.3 Proportion of droplets (named by HTO and genetic donor)

# proportion of genetically assigned droplets in each HTO
p3 <- ggcells(sce) + 
  geom_bar(
    aes(x = HTODemux.ID, fill = genetic_donor),
    position = position_fill(reverse = TRUE)) +
  coord_flip() +
  ylab("Frequency") +
  theme_cowplot(font_size = 8) + 
  scale_fill_manual(values = sce$colours$genetic_donor_colours)

# proportion of HTO assigned droplets in each genetic donor
p4 <- ggcells(sce) + 
  geom_bar(
    aes(x = genetic_donor, fill = HTODemux.ID),
    position = position_fill(reverse = TRUE)) +
  coord_flip() +
  ylab("Frequency") +
  theme_cowplot(font_size = 8) + 
  scale_fill_manual(values = sce$colours$hto_colours)

(p3 + p3 + facet_grid(~Capture) + plot_layout(widths = c(1, 2))) /
  (p4 + p4 + facet_grid(~Capture) + plot_layout(widths = c(1, 2))) +
  plot_layout(guides = "collect")

Version	Author	Date
29ed4f7	tcwong1994	2026-04-22

7.4 Number of droplets (named by sample ID)

# add colour
sampleID.HTO_colours <- setNames(
  c(palette.colors(nlevels(sce$sampleID.HTO)-2, "Tableau 10"), "#ced4da","#6c757d"),
  levels(sce$sampleID.HTO))
sce$colours$sampleID.HTO_colours <- sampleID.HTO_colours[sce$sampleID.HTO]

sampleID.genetics_colours <- setNames(
  c(palette.colors(nlevels(sce$sampleID.genetics)-2, "Set3"), "#d4a276","#9c6644"),
  levels(sce$sampleID.genetics))
sce$colours$sampleID.genetics_colours <- sampleID.genetics_colours[sce$sampleID.genetics]

# number of droplets assigned by HTO method
p5 <- ggcells(sce) + 
  geom_bar(aes(x = sampleID.HTO, fill = sampleID.HTO)) + 
  geom_text(stat='count', aes(x = sampleID.HTO, label=..count..), hjust=1, size=2) +
  coord_flip() + 
  ggtitle("By HTO") +
  ylab("Number of droplets") + 
  theme_cowplot(font_size = 8) + 
  scale_y_continuous(breaks=seq(0,16000,4000),limits=c(0,18000)) +
  scale_fill_manual(values = sce$colours$sampleID.HTO_colours) +
  theme(axis.title.y = element_blank(),
        axis.text.x=element_text(angle=45, hjust=1)) +
  guides(fill=guide_legend(title="sampleID"))

p5.facet <- ggcells(sce) + 
  geom_bar(aes(x = sampleID.HTO, fill = sampleID.HTO)) + 
  geom_text(stat='count', aes(x = sampleID.HTO, label=..count..), 
            hjust=1, size=2) +
  #ggtitle("By HTO") +
  ylab("Number of droplets") + 
  facet_grid(~Capture, scales = "fixed", space = "fixed") + 
  scale_y_continuous(breaks=seq(0,8000,2000), limits=c(0,10000)) +
  theme_cowplot(font_size = 8) + 
  scale_fill_manual(values = sce$colours$sampleID.HTO_colours) +
  theme(axis.title.y = element_blank(),
        axis.text.x=element_text(angle=45, hjust=1)) +
  guides(fill=guide_legend(title="sampleID")) +
    coord_flip()

# number of droplets assigned by genetic method
p6 <- ggcells(sce) + 
  geom_bar(aes(x = sampleID.genetics, fill = sampleID.genetics)) + 
  geom_text(stat='count', aes(x = sampleID.genetics, label=..count..), hjust=1, size=2) +
  coord_flip() + 
  ggtitle("By genetics") +
  ylab("Number of droplets") + 
  theme_cowplot(font_size = 8) + 
  scale_y_continuous(breaks=seq(0,16000,4000), limits = c(0,18000)) +
  scale_x_discrete(labels=gsub(".t[1-2]","",
                               levels(fct_drop(sce$sampleID.genetics)))) +
  scale_fill_manual(values = sce$colours$sampleID.genetics_colours) +
  theme(axis.title.y = element_blank(),
        axis.text.x=element_text(angle=45, hjust=1)) +
  guides(fill=FALSE)

p6.facet <- ggcells(sce) + 
  geom_bar(aes(x = sampleID.genetics, fill = sampleID.genetics)) + 
  geom_text(stat='count', aes(x = sampleID.genetics, label=..count..), 
            hjust=1, size=2) +
  ylab("Number of droplets") + 
  facet_grid(~Capture, space = "fixed", scales = "fixed") + 
  scale_y_continuous(breaks=seq(0,8000,2000), limits=c(0,10000)) +
  theme_cowplot(font_size = 8) + 
  scale_x_discrete(labels=gsub(".t[1-2]","",
                               levels(fct_drop(sce$sampleID.genetics)))) +
  scale_fill_manual(values = sce$colours$sampleID.genetics_colours) +
  theme(axis.title.y = element_blank(),
        axis.text.x=element_text(angle=45, hjust=1)) +
  guides(fill=guide_legend(title="sampleID")) +
  coord_flip()

(p5+p5.facet+plot_layout(width=c(1,2))) /
  (p6+p6.facet+plot_layout(width=c(1,2))) +
  plot_layout(guides="collect")

Version	Author	Date
29ed4f7	tcwong1994	2026-04-22

7.5 Proportion of droplets (named by sample ID)

# proportion of genetically assigned droplets in each HTO
p7 <- ggcells(sce) + 
  geom_bar(
    aes(x = sampleID.HTO, fill = sampleID.genetics),
    position = position_fill(reverse = TRUE)) +
  coord_flip() +
  ylab("Frequency") +
  theme_cowplot(font_size = 8) + 
  scale_fill_manual(values = sce$colours$sampleID.genetics_colours)

# proportion of HTO assigned droplets in each genetic donor
p8 <- ggcells(sce) + 
  geom_bar(
    aes(x = sampleID.genetics, fill = sampleID.HTO),
    position = position_fill(reverse = TRUE)) +
  coord_flip() +
  ylab("Frequency") +
  theme_cowplot(font_size = 8) + 
  scale_fill_manual(values = sce$colours$sampleID.HTO_colours)

((p7 + ggtitle("By HTO")) +
    p7 + facet_grid(~Capture) + plot_layout(widths = c(1, 2))) /
  ((p8 + ggtitle("By genetics")) + 
     p8 + facet_grid(~Capture) + plot_layout(widths = c(1, 2))) +
  plot_layout(guides = "collect")

Version	Author	Date
29ed4f7	tcwong1994	2026-04-22

8 Save object

demux_dir <- here("data","SCEs","demux")
if(!dir.exists(demux_dir)) {
  dir.create(demux_dir, recursive = TRUE)
}

out <- paste0(demux_dir,'/',
              paste0(batch_name,".cellbender.demux.SCE.rds"))

if(!file.exists(out)) saveRDS(sce, out)

sessionInfo()

R version 4.1.2 (2021-11-01)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: CentOS Linux 7 (Core)

Matrix products: default
BLAS:   /hpc/software/installed/R/4.1.2/lib64/R/lib/libRblas.so
LAPACK: /hpc/software/installed/R/4.1.2/lib64/R/lib/libRlapack.so

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       

attached base packages:
[1] stats4    stats     graphics  grDevices utils     datasets  methods  
[8] base     

other attached packages:
 [1] forcats_1.0.1               stringr_1.5.1              
 [3] janitor_2.1.0               vcfR_1.12.0                
 [5] dplyr_1.1.4                 scater_1.22.0              
 [7] scuttle_1.4.0               patchwork_1.2.0            
 [9] cowplot_1.1.3               SeuratObject_5.0.1         
[11] Seurat_4.4.0                ggplot2_3.5.2              
[13] here_1.0.1                  DropletUtils_1.14.2        
[15] SingleCellExperiment_1.16.0 SummarizedExperiment_1.24.0
[17] Biobase_2.54.0              GenomicRanges_1.46.1       
[19] GenomeInfoDb_1.30.1         IRanges_2.28.0             
[21] S4Vectors_0.32.4            BiocGenerics_0.40.0        
[23] MatrixGenerics_1.6.0        matrixStats_1.1.0          
[25] workflowr_1.7.0            

loaded via a namespace (and not attached):
  [1] utf8_1.2.4                spatstat.explore_3.2-7   
  [3] reticulate_1.36.1.9000    R.utils_2.12.2           
  [5] tidyselect_1.2.1          htmlwidgets_1.6.4        
  [7] grid_4.1.2                BiocParallel_1.28.3      
  [9] Rtsne_0.17                munsell_0.5.1            
 [11] ScaledMatrix_1.2.0        codetools_0.2-20         
 [13] ica_1.0-3                 statmod_1.5.0            
 [15] future_1.33.2             miniUI_0.1.1.1           
 [17] withr_3.0.0               spatstat.random_3.2-3    
 [19] colorspace_2.1-0          progressr_0.14.0         
 [21] highr_0.10                knitr_1.46               
 [23] rstudioapi_0.13           ROCR_1.0-11              
 [25] tensor_1.5.1              listenv_0.9.1            
 [27] labeling_0.4.3            git2r_0.31.0             
 [29] GenomeInfoDbData_1.2.11   polyclip_1.10-6          
 [31] pheatmap_1.0.12           farver_2.1.1             
 [33] rhdf5_2.38.1              rprojroot_2.1.1          
 [35] parallelly_1.37.1         vctrs_0.6.5              
 [37] generics_0.1.3            xfun_0.43                
 [39] R6_2.6.1                  ggbeeswarm_0.6.0         
 [41] rsvd_1.0.5                locfit_1.5-9.4           
 [43] bitops_1.0-9              rhdf5filters_1.6.0       
 [45] spatstat.utils_3.1-3      DelayedArray_0.20.0      
 [47] promises_1.3.0            scales_1.3.0             
 [49] pinfsc50_1.2.0            beeswarm_0.4.0           
 [51] gtable_0.3.5              beachmat_2.10.0          
 [53] globals_0.16.3            processx_3.8.0           
 [55] goftest_1.2-3             xaringanExtra_0.7.0      
 [57] spam_2.10-0               rlang_1.1.3              
 [59] splines_4.1.2             lazyeval_0.2.2           
 [61] spatstat.geom_3.2-9       yaml_2.3.8               
 [63] reshape2_1.4.4            abind_1.4-8              
 [65] httpuv_1.6.15             tools_4.1.2              
 [67] jquerylib_0.1.4           RColorBrewer_1.1-3       
 [69] ggridges_0.5.6            Rcpp_1.0.12              
 [71] plyr_1.8.9                sparseMatrixStats_1.6.0  
 [73] zlibbioc_1.48.2           purrr_1.2.1              
 [75] RCurl_1.98-1.9            ps_1.7.2                 
 [77] deldir_2.0-4              viridis_0.6.2            
 [79] pbapply_1.7-2             zoo_1.8-12               
 [81] ggrepel_0.9.1             cluster_2.1.8.1          
 [83] fs_1.6.4                  magrittr_2.0.3           
 [85] data.table_1.15.4         scattermore_1.2          
 [87] lmtest_0.9-40             RANN_2.6.1               
 [89] whisker_0.4               fitdistrplus_1.1-11      
 [91] mime_0.12                 evaluate_0.23            
 [93] xtable_1.8-4              gridExtra_2.3            
 [95] compiler_4.1.2            tibble_3.2.1             
 [97] KernSmooth_2.23-20        R.oo_1.24.0              
 [99] htmltools_0.5.8.1         mgcv_1.8-39              
[101] later_1.3.2               tidyr_1.3.1              
[103] lubridate_1.8.0           MASS_7.3-55              
[105] Matrix_1.6-5              permute_0.9-7            
[107] cli_3.6.2                 R.methodsS3_1.8.1        
[109] parallel_4.1.2            dotCall64_1.1-1          
[111] igraph_2.0.3              pkgconfig_2.0.3          
[113] getPass_0.2-2             sp_2.1-3                 
[115] plotly_4.10.4.9000        spatstat.sparse_3.0-3    
[117] memuse_4.2-1              vipor_0.4.5              
[119] bslib_0.3.1               dqrng_0.3.2              
[121] XVector_0.34.0            snakecase_0.11.0         
[123] callr_3.7.3               digest_0.6.35            
[125] sctransform_0.4.1         RcppAnnoy_0.0.22         
[127] vegan_2.5-7               spatstat.data_3.0-4      
[129] rmarkdown_2.26            leiden_0.4.3.1           
[131] uwot_0.1.14               edgeR_4.4.2              
[133] DelayedMatrixStats_1.16.0 shiny_1.8.1.1            
[135] lifecycle_1.0.4           nlme_3.1-155             
[137] jsonlite_1.8.8            Rhdf5lib_1.16.0          
[139] BiocNeighbors_1.12.0      viridisLite_0.4.2        
[141] limma_3.62.2              fansi_1.0.6              
[143] pillar_1.9.0              lattice_0.22-7           
[145] ggrastr_1.0.1             fastmap_1.1.1            
[147] httr_1.4.7                survival_3.2-13          
[149] glue_1.7.0                png_0.1-8                
[151] stringi_1.8.3             sass_0.4.9               
[153] HDF5Array_1.22.1          BiocSingular_1.10.0      
[155] ape_5.6-1                 irlba_2.3.5.1            
[157] future.apply_1.11.2

Demultiplexing the BAL pool: Pi1

Anson Wong

2026-04-22

1 Load packages

2 Set names for this pool

3 Read CellBender outputs

4 HTO demultiplexing

4.1 Prepare HTO data

4.2 HTO Demultiplexing

4.3 Add HTO demultiplexing results to SCE

5 Visualize HTO demultiplexing results

5.1 Number of singlets

Pi1-c1

Pi1-c2

5.2 HTO expression in ridge plot

Pi1-c1

Pi1-c2

5.3 HTO expression in heatmap

Pi1-c1

Pi1-c2

5.4 tSNE embeddings for HTOs

Pi1-c1

Pi1-c2

5.5 Number of demultiplexed droplets

Named by HTO

Named by sample ID

5.6 Proportion of singlets per capture

6 Genetic demultiplexing

6.1 Read vireo results

6.2 Summarise vireo results

6.3 Identify the best match between donors from the combination of captures

6.4 Add genetic demultiplexing results to SCE

7 Visualize genetic demultiplexing results

7.1 Table of number of droplets

7.2 Number of droplets (named by HTO and genetic donor)

7.3 Proportion of droplets (named by HTO and genetic donor)

7.4 Number of droplets (named by sample ID)

7.5 Proportion of droplets (named by sample ID)

8 Save object