Using ERCC spike-in features

Let us recall that the same set of ERCC spike-in features was added in identical quantity to each sample. Consequently, ERCC spike-in features should exhibit no biological variability, i.e., any variance in their counts should be technical in origin.

Let us fit a trend to the variance of spike-in features using the scran trendVar function; while accounting for the minor effect of technical factors described earlier. We do not block on the key experimental factors Time, Infection, and Status, as doing so may regress out biological effects of interest.

dm <- model.matrix(~Plate, data = colData(sce.norm))
colnames(dm)

## [1] "(Intercept)" "PlateP2"     "PlateP3"     "PlateP4"

Having defined the design matrix, we fit the mean-dependent trend to the gene-specific variances in the data set:

var.fit.ercc <- trendVar(
  sce.norm,  method = "loess", assay.type="logcounts", parametric = TRUE,
  design = dm, loess.args=list(span = 0.4), use.spikes=TRUE
)
names(var.fit.ercc)

## [1] "vars"     "means"    "resid.df" "design"   "trend"    "df2"

We may then decompose the gene-specific variance into biological and technical components using the ERCC spike-in trend as an estimate of technical variance:

var.out.ercc <- decomposeVar(
  sce.norm, var.fit.ercc, assay.type="logcounts", get.spikes=FALSE)
names(var.out.ercc)

## [1] "mean"    "total"   "bio"     "tech"    "p.value" "FDR"

Finally, we may plot the fitted trend, while highlighting ERCC spike-in features used to obtain it:

plot(var.out.ercc$mean, var.out.ercc$total, pch = 16, cex = 0.6,
  xlab = "Mean log-expression", ylab = "Variance of log-expression"
)
points(var.fit.ercc$mean, var.fit.ercc$var, col="red", pch=16)
o <- order(var.out.ercc$mean);
lines(var.out.ercc$mean[o], var.out.ercc$tech[o], col="red", lwd=2)

While the above figure suggests that the trend fits accurately spike-in features at extreme values the detection range, it also indicate substantial scatter around the trend for feature at intermediate levels of detection, and likely understimates technical variance in both ERCC spike-in and endogenous features at moderately high levels of detection (i.e., mean log-expression in the range 5-10).

Using endogenous features

Similarly to the use of ERCC spike-in features, let us use the scran trendVar method to estimate technical variation under the assumption stated above, accouting for the same technical factors of the experimental design:

var.fit.endo <- trendVar(
  sce.norm, method = "loess", parametric = TRUE,
  design = dm, loess.args=list(span = 0.1), assay.type="logcounts", use.spikes=FALSE
)
names(var.fit.endo)

## [1] "vars"     "means"    "resid.df" "design"   "trend"    "df2"

Similarly, the estimated variance is decomposed into its technical and biological components (using only the endogenous features processed above):

var.out.endo <- decomposeVar(
  sce.norm, var.fit.endo, assay.type="logcounts", get.spikes=FALSE)
names(var.out.endo)

## [1] "mean"    "total"   "bio"     "tech"    "p.value" "FDR"

In the figure below, we assess the suitability of the trend fitted to the endogenous variances by examining whether it is consistent with the spike-in variances.

plot(var.out.endo$mean, var.out.endo$total, pch=16, cex=0.6,
  xlab="Mean log-expression", ylab="Variance of log-expression"
)
o <- order(var.out.endo$mean)
lines(var.out.endo$mean[o], var.out.endo$tech[o], col="dodgerblue", lwd=2)
spike.fit <- trendVar(sce.norm, assay.type="logcounts", use.spikes=TRUE) # To compute spike-in variances.
points(spike.fit$mean, spike.fit$var, col="red", pch=16)

Notably, the trend passes through or close to most of the spike-in variances, suggesting that the assumption (that most genes have low levels of biological variability) is suitable for this data set. Notably, this strategy exploits the large number of endogenous genes to obtain a stable trend, with the ERCC spike-in features primarily used as diagnostic features rather than in the trend fitting itself.

Importantly, if the assumption (that most genes have low levels of biological variability) did not hold, one would rather use the trend fitted to the ERCC spike-in features using the default use.spikes=TRUE settings. Notably, this would sacrifice stability to reduce systematic errors in the estimate of the biological component for each gene (Lun, McCarthy, and Marioni 2016).

Identification

hvg.out <- var.out.endo[with(var.out.endo, which(FDR <= 0.05 & bio >= 0.5)),]
hvg.out <- hvg.out[order(hvg.out$bio, decreasing=TRUE),]

The above criteria have identified 1270 endogenous genes as highly variable.

The 100 genes with highest biological component of variance may be examined in the second tab.

topHvgIds <- rownames(hvg.out)[1:10]

Time

The figure above generally indicates a progressive up-regulation of gene expression over time.

Infection

## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.

The figure above further links the high variance of those genes to the differences between bacteria-stimulated and mock-infected cells, particularly at the later time-points.

Status

## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.

The figure above reinforces clear-cut differences between cells stimulated by bacterial infection or exposure and uninfected cells, particularly at the later time-points.

Computation

The initial step of the goseq procedure is to create a named vector that indicates whether each gene (i.e., names) is highly variable (1) or not (0):

endoFeatures <- rownames(sce.norm)[!isSpike(sce.norm, "ERCC")]
hvgs <- as.integer(endoFeatures %in% rownames(hvg.out))
names(hvgs) <- endoFeatures
table(hvgs)

## hvgs
##    0    1 
## 6956 1270

Next, goseq is capable of accounting for gene length bias during GO testing. In the current case, the gene length information reported by featureCounts is used:

geneLengths <- width(rowRanges(sce.norm[rownames(hvg.out),]))
names(geneLengths) <- rownames(hvg.out)

Using the above gene length data, the probability weighting function (PWF) is calculated for the full gene set in this analysis:

pwf <- goseq::nullp(
  hvgs, bias.data = width(rowRanges(sce.norm[names(hvgs),])), plot.fit=TRUE)

Finally, the goseq function is used to test for enrichment amongst the HVGs:

GO.out <- goseq::goseq(pwf, "hg38", "ensGene")

## Fetching GO annotations...

## Loading required package: AnnotationDbi

## 

## For 744 genes, we could not find any categories. These genes will be excluded.

## To force their use, please run with use_genes_without_cat=TRUE (see documentation).

## This was the default behavior for version 1.15.1 and earlier.

## Calculating the p-values...

## 'select()' returned 1:1 mapping between keys and columns

The ontology field (i.e., "BP", "MF", or "CC") is reformatted as a factor:

GO.out$ontology <- as.factor(GO.out$ontology)

Euclidian distance

Estimate distance between cells using Euclidian distance:

d.e <- dist(t(norm.endo.hvg), diag = TRUE, upper = TRUE)
mat.e <- as.matrix(d.e)

We may then cluster samples according to their Euclidian distance, using Ward’s criterion to minimize the total variance within each cluster (Lun, McCarthy, and Marioni 2016):

h.e <- hclust(d.e, method = "ward.D2")
ord.e <- mat.e[h.e$order, h.e$order]
rm(mat.e)

Assemble the heat map:

ht_column <- HeatmapAnnotation(
  df = colData(sce.endo)[h.e$order,c("Time", "Status", "Infection")],
  col = list(
    Time = col.time,
    Status = col.status,
    Infection = col.infection
  )
)
ht_e <- Heatmap(
  ord.e, name = "distance", column_title = "Euclidian distance",
  top_annotation = ht_column,
  cluster_rows = FALSE, cluster_columns = FALSE,
  show_row_names = FALSE, show_column_names = FALSE
)

Spearman correlation

Estimate distance between cells using Spearman’s \(\rho\) correlation.

c.s <- MKmisc::corDist(t(norm.endo.hvg), "spearman", diag = TRUE, upper = TRUE)
mat.c.s <- as.matrix(c.s)

Cluster samples:

h.c.s <- hclust(c.s)
ord.c.s <- mat.c.s[h.c.s$order, h.c.s$order]
rm(c.s)

Assemble the heat map:

ht_column <- HeatmapAnnotation(
  df = colData(sce.endo)[h.c.s$order,c("Time", "Status", "Infection")],
  col = list(
    Time = col.time,
    Status = col.status,
    Infection = col.infection
  )
)
ht_s <- Heatmap(
  ord.c.s, name = "distance", column_title = "Spearman absolute correlation distance",
  top_annotation = ht_column,
  cluster_rows = FALSE, cluster_columns = FALSE,
  show_row_names = FALSE, show_column_names = FALSE
)

Pearson correlation

Estimate distance between cells using Pearson’s \(r\) correlation:

c.p <- MKmisc::corDist(t(norm.endo.hvg), diag = TRUE, upper = TRUE)
mat.c.p <- as.matrix(c.p)

Cluster samples using Pearson’s correlation:

h.c.p <- hclust(c.p)
ord.c.p <- mat.c.p[h.c.p$order, h.c.p$order]
rm(c.p)

Assemble the heat map:

hn_column <- HeatmapAnnotation(
  df = colData(sce.endo)[h.c.p$order,c("Time", "Status", "Infection")],
  col = list(
    Time = col.time,
    Status = col.status,
    Infection = col.infection
  )
)
ht_p <- Heatmap(
  ord.c.p, name = "distance", column_title = "Pearson absolute correlation distance",
  top_annotation = hn_column,
  cluster_rows = FALSE, cluster_columns = FALSE,
  show_row_names = FALSE, show_column_names = FALSE
)

Euclidian

Spearman

Pearson