Open xie186 opened 3 years ago
HelenaLC
commented
3 years ago Thank you for the kind words!
You said you checked the sub object, but your line above has sce in it? Could you maybe print bout sce and sub here so I can have a look what these looks like?
Yes, compexprs contain compensated normalized expression values - assuming that you ran compCytof(..., assay = "normcounts") (i.e., running the compensation on the normalized data, not the default assay = "counts").
xie186
commented
3 years ago Thank you for your response. This is the print output:
> print(sce)
class: SingleCellExperiment
dim: 63 337525
metadata(4): experiment_info bc_key sep_cutoffs mhl_cutoff
assays(7): counts exprs ... compcounts compexprs
rownames(63): 75As CD15 ... 208Pb CD45
rowData names(6): channel_name marker_name ... bead_ch is_bc
colnames: NULL
colData names(6): file_id remove ... delta mhl_dist
reducedDimNames(0):
altExpNames(0):
> print(sub)
class: SingleCellExperiment
dim: 3 337525
metadata(4): experiment_info bc_key sep_cutoffs mhl_cutoff
assays(7): counts exprs ... compcounts compexprs
rownames(3): DNA1 DNA2 dead
rowData names(6): channel_name marker_name ... bead_ch is_bc
colnames: NULL
colData names(2): bc_id i
reducedDimNames(0):
altExpNames(0):I also pasted here my code. Basically the code is from your paper:
#' ## Load library
#'
library(CATALYST)
library(dplyr)
library(flowCore)
library(flowWorkspace)
library(mvtnorm)
library(openCyto)
library(ggcyto)
library(ggplot2)
qc_theme <- list(
theme_bw(base_size = 8), theme(
panel.grid.minor = element_blank(),
panel.grid.major.x = element_blank(),
plot.title = element_text(face = "bold"),
axis.text = element_text(color = "black"),
axis.text.x = element_text(angle = 45, hjust = 1)))
#' ## Read `FCS` files
#'
#' A SCE can be constructed using CATALYST’s `prepData()` function, which accepts
#' a path to a directory with one or many FCS files, a character vector of FCS
#' filenames, a single or list of flowFrame(s), or a
#' https://bioconductor.org/packages/3.11/bioc/html/flowCore.html (flowCore package24).
#' By default (transform = TRUE), an arcsinh-transformation with cofactor = 5 is applied to
#' the input (count) data, and the resulting expression matrix is stored in the exprs assay
#' slot of the output SCE:
# construct ’SingleCellExperiment’
data_path = "/depot/bioinf/data/Projects/CyTOF/test_data/"
setwd("/depot/bioinf/data/Projects/CyTOF/test_data/")
fcs <- list.files(data_path, "acquisition", full.names = TRUE)
fcs
(sce <- prepData(fcs, transform = TRUE, cofactor = 5))
## The cofactor used for transformation
(int_metadata(sce)$cofactor)
## dimmencial of counts
class(sce@assays@data$counts)
dim(sce@assays@data$counts)
## cofactor-5 arcsinh-transformed counts
class(sce@assays@data$exprs)
dim(sce@assays@data$exprs)
## How many cells in each file
head(colData(sce))
table(colData(sce))
#cell_num <- data.frame(
# file_id = levels(sce$file_id),
# n_cells = tabulate(sce$file_id))
#cell_num
## we rename the cell metadata variable sample_id to file_id to avoid ambiguity:
i <- match("sample_id", names(colData(sce)))
names(colData(sce))[i] <- "file_id"
# store character vector or channels names
chs <- channels(sce)
# store DNA & live channel indices
dna <- grep("^Ir", chs)
live <- grep("191|194", chs)
#' ## Filtering for stable signal
# plot channels of interest vs. time
coi <- chs[c(dna[1], which(rowData(sce)$use_channel))]
coi
plotScatter(sce, chs = c("Time", coi), label = "both") +
labs(y = "expression") +
scale_x_continuous(
expression("Time ("*10^6~"ms)"),
labels = function(u) u/1e6) +
theme_bw(base_size = 8) + theme(
aspect.ratio = 2/3,
panel.grid = element_blank(),
axis.text = element_text(color = "black"),
strip.background = element_rect(fill = NA))
#' ## Data filtering
## Let's assume we want to remove 0 to 1 10^6 ms
# construct ’GatingSet’
ff <- sce2fcs(sce[dna, ], assay = "exprs")
gs <- GatingSet(flowSet(ff))
# apply rectangular gate to exclude unstable signal
min_t <- 0
max_t <- 1000000
gs_add_gating_method(
gs, alias = "stable",
pop = "-", parent = "root",
dims = paste0("Time,", chs[dna[1]]),
gating_method = "boundary",
gating_args = sprintf("min=c(%s,0),max=c(%s,10)", min_t, max_t))
# plot scatter of DNA vs. Time
ggcyto(gs,
aes_string("Time", chs[dna[1]])) +
geom_hex(bins = 128) +
geom_gate("stable") +
facet_null() + theme_bw() +
ggtitle(NULL) + theme(
legend.position = "none",
panel.grid = element_blank())
# subset selected events
sce_fil <- sce[, gh_pop_get_indices(gs, "stable")]
plotScatter(sce_fil, chs = c("Time", coi), label = "both") +
labs(y = "expression") +
scale_x_continuous(
expression("Time ("*10^6~"ms)"),
labels = function(u) u/1e6) +
theme_bw(base_size = 8) + theme(
aspect.ratio = 2/3,
panel.grid = element_blank(),
axis.text = element_text(color = "black"),
strip.background = element_rect(fill = NA))
#' ## Normalization
# specify path to reference beads
ref_beads <- file.path(data_path, "normalization_beads.fcs")
# apply bead-based normalization
system.time(res <- normCytof(sce,
beads = "dvs", norm_to = ref_beads,
remove_beads = TRUE, overwrite = FALSE))
#' When remove_beads = TRUE (the default), normCytof() will return a list of three SCEs containing filtered,
#' `bead` and `remove` events, respectively, as well as two ggplot objects:
names(res)
table(res$data$is_bead)
table(res$data$remove)
table(res$data$file_id)
# view no. of remaining, bead & removed events
ns <- sapply(res[1:3], ncol)
ps <- sprintf("%1.2f", ns/ncol(sce)*100)
stat_cell_num <- data.frame(t(cbind("# events" = ns, "% of total" = ps)))
stat_cell_num
#' As a first quality control plot, `res$scatter` renders scatter plots of bead channels (x-axis)
#' versus DNA (y-axis), where events identified as beads as well as their expression range are
#' highlighted in color; bead events should have low DNA intensity (since they are not cells) and
#' high intensities across all bead channels.
res$scatter
res$lines
# compute mean bead channel counts for current run
(rowData(res$beads))
is_bead <- rowData(res$beads)$bead_ch # get bead channels
bead_cs <- counts(res$beads)[is_bead, ] # subset counts
rownames(bead_cs) <- chs[is_bead] # use channels as names
(bead_ms <- rowMeans(bead_cs)) # compute means
# read in reference mean bead channel counts
ref <- read.csv(file.path(data_path, "ref_bead_counts.csv"))
# join into single tidy data.frame
df <- bind_rows(ref, bead_ms, .id = "group")
library(reshape2)
df <- melt(df, id.var = "group")
# boxplot of reference vs. current run’s mean bead channel counts
ggplot(df, aes(variable, value)) +
geom_boxplot(data = df[df$group == 1, ]) +
geom_point(data = df[df$group == 2, ],
col = "red", pch = 4, stroke = 1) +
labs(x = "bead channel", y = "mean count") +
qc_theme + ggtitle(
"QC on bead channel counts",
"[-] = reference | x = current run")
#' After normalization, we overwrite the input dataset with the filtered subset
#' that no longer includes bead events, or bead-bead and bead-cell doublets:
sce <- res$data
#' ## Debarcoding
# read in debarcoding scheme
fn <- file.path(data_path, "debarcoding_scheme.csv")
bc_key <- read.csv(fn, row.names = 1, check.names = FALSE)
# all barcodes are positive for exactly 3 barcoding channels
all(rowSums(bc_key) == 3)
# remove empty barcodes from debarcoding scheme
is_empty <- grepl("empty", rownames(bc_key))
bc_key <- bc_key[!is_empty, ]
bc_ids <- rownames(bc_key)
# do preliminary barcode assignments
system.time(sce <- assignPrelim(sce, bc_key, assay = "normexprs"))
# view barcode channels
channels(sce)[rowData(sce)$is_bc]
table(rowData(sce)$is_bc)
# tabulate number of (un)assigned events
tab_bc_id = table(sce$bc_id == 0)
tab_bc_id[2]/(tab_bc_id[1]+tab_bc_id[2])
#' ## Estimation of separation cutoffs
plotYields(sce, which = "0")
sce <- estCutoffs(sce)
metadata(sce)$sep_cutoffs
plotYields(sce, which = "PBMC_R1")
# store preliminary barcode IDs
bc_ids0 <- sce$bc_id
# apply global & sample-specific separation cutoff(s)
sce_glob <- applyCutoffs(sce, sep_cutoffs = 0.15, mhl_cutoff = 30)
sce_spec <- applyCutoffs(sce, mhl_cutoff = 30)
# compare cell yields for both cutoff strategies
c(global = mean(sce_glob$bc_id == 0),
specific = mean(sce_spec$bc_id == 0))
# proceed with sample-specific filtering
sce <- sce_spec
# compute number of events per population
# before vs. after applying separation cutoffs
barplot(rbind(table(bc_ids0), table(sce$bc_id)),
beside = TRUE, ylab = "cell count",
las = 2, cex.axis = 0.5, cex.names = 0.5)
legend("topright", fill = c("black", "grey"),
legend = c("before filtering", "after filtering"))
# event plots for unassigned events
# & barcode population D1
plotEvents(sce, which = c(0, "Tumor_S1", "Tumor_S2"), n = 25)
# plotMahal(sce, which = "0")
#' ## Compensation
# read in pre-computed spillover matrix
sm <- file.path(data_path, "spillover_matrix.csv")
sm <- read.csv(sm, row.names = 1)
# apply NNLS compensation
system.time(
sce <- compCytof(sce, sm, method = "nnls",
assay = "normcounts", overwrite = FALSE))
#' To visually inspect how compensation affects signal intensities, we can generate scatter plots
#' pre- and post-compensation; an exemplary pair of channels is shown in the figure. In such a plot,
#' we can observe a slight positive association between the signals of spill-affected channels,
#' which should be removed upon compensation.
i <- grep("173|174", chs, value = TRUE)
p1 <- plotScatter(sce,
chs = i,
label = "channel",
assay = "normexprs") +
ggtitle("Uncompensated")
p2 <- plotScatter(sce,
chs = i,
label = "channel",
assay = "compexprs") +
ggtitle("Compensated") +
ylab(NULL)
library(patchwork)
wrap_plots(p1, p2)
#' ## Gating
#'
.scatter <- function(gs, chs, gate_id = NULL,
subset = ifelse(is.null(gate_id), "root", "_parent_")) {
p <- ggcyto(gs, max_nrow_to_plot = 1e5,
aes_string(chs[1], chs[2]), subset) +
geom_hex(bins = 100) + facet_wrap(~ name, ncol = 5) +
(if (is.null(gate_id)) list() else geom_gate(gate_id)) +
ggtitle(NULL) + theme_bw(base_size = 8) + theme(
aspect.ratio = 1,
legend.position = "none",
panel.grid.minor = element_blank(),
strip.background = element_rect(fill = NA),
axis.text = element_text(color = "black"),
axis.text.x = element_text(angle = 45, hjust = 1))
suppressMessages(p + coord_equal(expand = FALSE,
xlim = c(-1, 11), ylim = c(-1, 11)))
}
#' ### Gating on cells
# subset DNA & live channels
sub <- sce[union(dna, live), ]
# add metadata variable ’i’ to track cell indices
colData(sub) <- DataFrame(
bc_id = sub$bc_id,
i = seq_len(ncol(sce)))
## split SCE by sample
fs <- sce2fcs(sub,
assay = "compexprs",
split_by = "bc_id",
keep_cd = TRUE)
Thank you for the nice work flow published on F1000: https://f1000research.com/articles/9-1263/v1
Question 1:
I am able to successfully run the R code upto the code below:
The code above give me an error: Error: Argument 'exprs' must be numeric matrix with colnames attribute set
I checked
subobject.The other two judgement were all
TRUE. What's your suggestion toQuestion 2: My understanding is the values in
compexprsis compensated and normalizedexprand if I want to do clustering and tSNE/UMAP, I should use the valuescompexprsafter gating. Please let me know whether my understanding is correct? Thanks.Here is my session info: