Preprocessing

Let us first calculate the log₂-transformed count per million (CPM) after adding a default \(0.25\) count to all genes, to avoid computing the logarithm of \(0\):

logCPM.cells <- edgeR::cpm(counts(sce.pass), log=TRUE, prior.count=0.25)
logCPM.bulk <- edgeR::cpm(counts(sce.bulk), log=TRUE, prior.count=0.25)

Let us progressively collect data for each experimental group:

bulkSCdata <- data.frame()
for (x_group in levels(sce.pass$Group)){
  message(x_group)
  time <- gsub("^(.*h)_.*", "\\1", x_group)
  infection <- gsub(".*_(.*)_.*", "\\1", x_group)
  status <- gsub(".*_.*_(.*)", "\\1", x_group)
  treatment <- gsub(".*_(.*_.*)", "\\1", x_group)
  cellIndex <- (sce.pass$Group == x_group)
  featureIndex <- (
    rowSums(counts(sce.pass)) > 0 |
      rowSums(counts(sce.bulk)) > 0)
  cellAvg <- rowMeans(logCPM.cells[featureIndex, cellIndex])
  bulkExpr <- logCPM.bulk[
    featureIndex,
    which(sce.bulk$Time == time & sce.bulk$Infection == infection)
  ]
  geneNames <- names(cellAvg)
  bulkSCdata <- rbind(
    bulkSCdata,
    data.frame(
      Cells = cellAvg,
      Bulk = bulkExpr[geneNames], # ensure order identical to single cells
      Gene = geneNames,
      Group = x_group, Treatment = treatment,
      Time = time, Infection = infection, Status = status
    )
  )
}

Stimulated bulk samples are expected to be a mixture of infected (minority; < 5%) and exposed cells (vast majority; > 95%). We may therefore compare the average expression in each subset of cells against the bulk average, both in log₂(CPM), to identify potential differences. For reference, let us also indicate the identity line (\(y=x\)) in grey:

bulkSCdata$ERCC <- grepl("^ERCC-", bulkSCdata$Gene)
ggplot(bulkSCdata) +
  geom_hex(aes(Bulk, Cells), subset(bulkSCdata, !ERCC)) +
  geom_point(aes(Bulk, Cells), subset(bulkSCdata, ERCC),size=0.1,colour="red")+
  geom_abline(slope = 1, intercept = 0, colour = "grey") +
  coord_fixed(ratio = 1) + labs(y="Cell average") + guides(fill = FALSE) +
  facet_grid(Infection + Status ~ Time) +
  theme_minimal() + theme(
    axis.text.x = element_text(size = rel(0.5)),
    axis.text.y = element_text(size = rel(0.5))
  ) + guides(fill = "none")

The figure above displays generally higher detection levels in bulk samples for a large proportion of genes; this demonstrates the higher frequency of dropout events in single cell samples, and the average measurement of bulk that rescues detection of rarer transcripts expressed stochastically at low levels or in few cells.

Counts per million (CPM)

First, let us initalise the matrix of CPM:

SingleCellExperiment::cpm(sce.pass) <- edgeR::cpm(counts(sce.pass))
range(colSums(cpm(sce.pass)))

## [1] 1e+06 1e+06

Note: The SingleCellExperiment::cpm method is an accessor to the corresponding slot of SingleCellExperiment objects, while the edgeR::cpm method calculates normalised counts per million (CPM) for given a matrix of read counts. Furthermore, in the absence of size factors, CPM are computed using the total count of read (i.e., library size) for each sample, without further normalisation of library sizes between samples.

Detection frequency and mean expression level

For this purpose, let us first calculate updated QC metrics–in the subset of cells that passed quality control in the previous section–using once more the scater calculateQCMetrics method:

sce.pass <- calculateQCMetrics(
  sce.pass, feature_controls = list(ERCC = ERCCs, MT = MTs)
)

## Note that the names of some metrics have changed, see 'Renamed metrics' in ?calculateQCMetrics.
## Old names are currently maintained for back-compatibility, but may be removed in future releases.

We may then examine for each feature the proportion of cells in which it is detected (non-zero expression level) against its average expression level across all cells:

plotExprsFreqVsMean(
  sce.pass,
  feature_controls = rowData(sce.pass)[,"is_feature_control_ERCC"]
)

## Warning in .get_all_sf_sets(object): spike-in set 'ERCC' should have its
## own size factors

## `geom_smooth()` using method = 'loess'

Spike-in expectations

Let us compare the expected number of spike-in molecules introduced in each sample (calculated in this earlier section) against the average count detected across all single cells, while indicating the identity relationship (\(y=x\)) as a solid black line:

mixC$avgCount <- rowMeans(counts(sce.pass)[mixC$ID,])
ggplot(mixC) +
  geom_point(aes(ERCCmolecules, avgCount)) +
  geom_abline(slope = 1, intercept = 0, linetype = "dashed", colour = "red") +
  labs(
    x = "Estimated ERCC molecules injected",
    y = "Average observed fragment count in cells passed QC") +
  scale_x_log10() + scale_y_log10() + coord_fixed(ratio = 1) + theme_minimal()

Count of cells with detectable expression

We may compare the above figure with the a smoother representation of the count of cells in which each feature is detected against the average log₁₀-transformed count of reads assigned to that same feature, while marking the same cut-off threshold of an average 1 count across all cells:

numcells <- nexprs(sce.pass, byrow = TRUE)
is.ercc <- rowData(sce.pass)[,"is_feature_control_ERCC"]
smoothScatter(
  log10(ave.cpm), numcells, xlab=expression(Log[10]~"average CPM"),
  ylab="Number of expressing cells"
)
points(
  log10(ave.cpm[is.ercc]), numcells[is.ercc],
  col = "red", pch = 16, cex = 0.5
)
abline(v = log10(1), col = "blue", lwd = 2, lty = 2)