5 Preprocessing
5.1 Control data quality
Remove variables (genes) and samples (cells) with low quality, by processing the following three steps:
- remove variables based on expression profiles across samples,
- remove samples based on the numbers of reads and nonzero expressed variables,
- remove variables based on the mean read counts across samples.
First of all, add metadata for both variables and samples using
ASURAT function add_metadata().
sclc <- add_metadata(sce = sclc, mitochondria_symbol = "^MT-")Examine the expression levels of known SCLC marker genes (Ireland, et al., 2020), namely ASCL1, NEUROD1, YAP1, and POU2F3.
genes <- c("ASCL1", "NEUROD1", "YAP1", "POU2F3")
sce <- sclc[, which(colData(sclc)$nReads > 2000)]
set.seed(1)
inds <- sample(ncol(sce), size = 1000, replace = FALSE)
subsce <- sce[genes, inds]
mat <- log(as.matrix(assay(subsce, "counts")) + 1)
set.seed(1)
filename <- "figures/figure_01_0005.png"
png(file = filename, height = 200, width = 520, res = 100)
#png(file = filename, height = 580, width = 1600, res = 300)
p <- ComplexHeatmap::Heatmap(mat, column_title = "SCLC",
name = "Log1p\nexpression", cluster_rows = FALSE,
show_row_names = TRUE, row_names_side = "right",
show_row_dend = FALSE, show_column_names = FALSE,
column_dend_side = "top", show_parent_dend_line = FALSE)
p
dev.off()
# mtx <- t(colData(subsce)$nReads) ; rownames(mtx) <- "nReads"
# q <- ComplexHeatmap::Heatmap(mtx, name = "nReads", show_row_names = TRUE,
# row_names_side = "right", show_row_dend = FALSE,
# show_column_names = FALSE, show_column_dend = FALSE,
# col = circlize::colorRamp2(c(min(mtx), max(mtx)),
# c("cyan", "magenta")))
# p <- p %v% q
# p
# dev.off()
5.1.1 Remove variables based on expression profiles
ASURAT function remove_variables() removes variable
(gene) data such that the numbers of non-zero expressing samples (cells)
are less than min_nsamples.
sclc <- remove_variables(sce = sclc, min_nsamples = 10)5.1.2 Remove samples based on expression profiles
Qualities of sample (cell) data are confirmed based on proper
visualization of colData(sce).
title <- "SCLC"
df <- data.frame(x = colData(sclc)$nReads, y = colData(sclc)$nGenes)
p <- ggplot2::ggplot() +
ggplot2::geom_point(ggplot2::aes(x = df$x, y = df$y), size = 1, alpha = 1) +
ggplot2::labs(title = title, x = "Number of reads", y = "Number of genes") +
ggplot2::theme_classic(base_size = 20) +
ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5, size = 20))
filename <- "figures/figure_01_0010.png"
ggplot2::ggsave(file = filename, plot = p, dpi = 50, width = 5, height = 5)
df <- data.frame(x = colData(sclc)$nReads, y = colData(sclc)$percMT)
p <- ggplot2::ggplot() +
ggplot2::geom_point(ggplot2::aes(x = df$x, y = df$y), size = 1, alpha = 1) +
ggplot2::labs(title = title, x = "Number of reads", y = "Perc of MT reads") +
ggplot2::theme_classic(base_size = 20) +
ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5, size = 20))
filename <- "figures/figure_01_0011.png"
ggplot2::ggsave(file = filename, plot = p, dpi = 50, width = 5, height = 5)
ASURAT function remove_samples() removes sample (cell)
data by setting cutoff values for the metadata.
sclc <- remove_samples(sce = sclc, min_nReads = 1400, max_nReads = 40000,
min_nGenes = 1150, max_nGenes = 1e+10,
min_percMT = 0, max_percMT = 15)5.1.3 Remove variables based on the mean read counts
Qualities of variable (gene) data are confirmed based on proper
visualization of rowData(sce).
title <- "SCLC"
aveexp <- apply(as.matrix(assay(sclc, "counts")), 1, mean)
df <- data.frame(x = seq_len(nrow(rowData(sclc))),
y = sort(aveexp, decreasing = TRUE))
p <- ggplot2::ggplot() + ggplot2::scale_y_log10() +
ggplot2::geom_point(ggplot2::aes(x = df$x, y = df$y), size = 1, alpha = 1) +
ggplot2::labs(title = title, x = "Rank of genes", y = "Mean read counts") +
ggplot2::theme_classic(base_size = 20) +
ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5, size = 20))
filename <- "figures/figure_01_0015.png"
ggplot2::ggsave(file = filename, plot = p, dpi = 50, width = 5, height = 5)
ASURAT function remove_variables_second() removes
variable (gene) data such that the mean read counts across samples are
less than min_meannReads.
sclc <- remove_variables_second(sce = sclc, min_meannReads = 0.20)dim(sclc)[1] 6346 22835.2 Normalize data
Perform bayNorm() (Tang et al., Bioinformatics, 2020)
for attenuating technical biases with respect to zero inflation and
variation of capture efficiencies between samples (cells).
mat <- as.matrix(assay(sclc, "counts"))
BETA <- bayNorm::BetaFun(Data = mat, MeanBETA = 0.06)
bayout <- bayNorm::bayNorm(mat, BETA_vec = BETA[["BETA"]], mode_version = TRUE)
assay(sclc, "normalized") <- bayout$Bay_outPerform log-normalization with a pseudo count.
assay(sclc, "logcounts") <- log(assay(sclc, "normalized") + 1)Center row data.
mat <- assay(sclc, "logcounts")
assay(sclc, "centered") <- sweep(mat, 1, apply(mat, 1, mean), FUN = "-")Set gene expression data into altExp(sce).
sname <- "logcounts"
altExp(sclc, sname) <- SummarizedExperiment(list(counts = assay(sclc, sname)))Add ENTREZ Gene IDs to rowData(sce).
dictionary <- AnnotationDbi::select(org.Hs.eg.db::org.Hs.eg.db,
key = rownames(sclc),
columns = "ENTREZID", keytype = "SYMBOL")
dictionary <- dictionary[!duplicated(dictionary$SYMBOL), ]
rowData(sclc)$geneID <- dictionary$ENTREZID
