4.1 Processing a subset of the HCA bone marrow data

For the scope of the vignette, we subset some cells in the bone marrow dataset to reduce the runtime, and apply some of the steps one would ideally follow in analysing droplet based datasets. For more information on how to properly process such datasets, please refer to the amazing set of resources available in the simpleSingleCell workflow package.

In brief: we start loading the required libraries for preprocessing, and taking a subset of the bone marrow dataset

library(scran)
library(scater)

set.seed(42)
sce <- sce_bonemarrow[, sample(seq_len(ncol(sce_bonemarrow)), 1000, replace = FALSE)]

First, we relabel the rows with the gene symbols for easier reading with the uniquifyFeatureNames() function from scater. Then we compute some QC metrics, using calculateQCMetrics().

rownames(sce) <- uniquifyFeatureNames(rowData(sce)$ID, rowData(sce)$Symbol)
head(rownames(sce))
#> [1] "RP11-34P13.3"  "FAM138A"       "OR4F5"         "RP11-34P13.7" 
#> [5] "RP11-34P13.8"  "RP11-34P13.14"

counts(sce) <- as.matrix(counts(sce))
sce <- scater::calculateQCMetrics(sce)

We proceed with normalization, performed with the deconvolution method implemented in scran - a pre-clustering step is done in advance, to avoid pooling cells that are very different between each other.

set.seed(42)
clusters <- quickCluster(sce, method="igraph", min.mean=0.1, irlba.args=list(maxit=1000))
table(clusters)
#> clusters
#>   1   2   3   4   5 
#> 222 210 129 178 261
sce <- computeSumFactors(sce, min.mean=0.1, cluster=clusters)
summary(sizeFactors(sce))
#>      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
#>  0.000465  0.325936  0.533787  1.000000  0.757043 28.645950
sce <- normalize(sce)

In the following lines of code, we model the mean-variance trend (with technical noise as Poisson), to extract the biological component of the variance for the genes under inspection.

new.trend <- makeTechTrend(x=sce)
fit <- trendVar(sce, use.spikes=FALSE, loess.args=list(span=0.05))
# plot(fit$mean, fit$var, pch=16)
# curve(fit$trend(x), col="dodgerblue", add=TRUE)
# curve(new.trend(x), col="red", add=TRUE)

fit0 <- fit
fit$trend <- new.trend
dec <- decomposeVar(fit=fit)
top.dec <- dec[order(dec$bio, decreasing=TRUE),] 
head(top.dec)
#> DataFrame with 6 rows and 6 columns
#>                    mean            total              bio              tech
#>               <numeric>        <numeric>        <numeric>         <numeric>
#> S100A8 1.11493057349193 5.38626448392199 4.29483312079354  1.09143136312844
#> S100A9 1.10437743825487 5.16160825818039 4.07615753701908  1.08545072116131
#> MALAT1  7.3308112726242 4.61018289683281 3.92709581513978 0.683087081693032
#> LYZ    1.15300320649838 4.99483962106588 3.88242885864087  1.11241076242501
#> HBB    1.08439490889895 4.24429550236879 3.17036827050413  1.07392723186466
#> RPS27  5.05339043152688 3.87882156515147 2.66980093676179  1.20902062838968
#>                      p.value                   FDR
#>                    <numeric>             <numeric>
#> S100A8                     0                     0
#> S100A9                     0                     0
#> MALAT1                     0                     0
#> LYZ                        0                     0
#> HBB                        0                     0
#> RPS27  5.68030787893909e-229 3.54430173468469e-226
plotExpression(sce, features=rownames(top.dec)[1:10])

With the denoisePCA function from scran, we can compute a reduced dimension view of the data, accounting for the Poisson technical trend.

Once we computed the PCA, we provide that as initialization to runTSNE, and obtain the t-SNE results, stored in the appropriate slot of our sce object.

set.seed(42)
sce <- denoisePCA(sce, technical=new.trend, approximate=TRUE)
ncol(reducedDim(sce, "PCA"))
#> [1] 7
plot(attr(reducedDim(sce), "percentVar"), xlab="PC",
     ylab="Proportion of variance explained")
abline(v=ncol(reducedDim(sce, "PCA")), lty=2, col="red")

plotPCA(sce, ncomponents=3, colour_by="log10_total_features_by_counts")

set.seed(42)
sce <- runTSNE(sce, use_dimred="PCA", perplexity=30)
plotTSNE(sce, colour_by="log10_total_features_by_counts")

After this, we can compute a shared nearest neighbour graph to identify clusters in our dataset. Once the cluster memberships are defined, we assign this vector to a corresponding colData slot, and plot the t-SNE for our subset, coloured accordingly.

snn.gr <- buildSNNGraph(sce, use.dimred="PCA")
clusters <- igraph::cluster_walktrap(snn.gr)
sce$Cluster <- factor(clusters$membership)
table(sce$Cluster)
#> 
#>   1   2   3   4   5   6   7   8   9  10 
#> 183 171  69 240  69  43 106  35  37  47
plotTSNE(sce, colour_by="Cluster")

Even if we only took a subset of the available full data, we can observe that different clusters are indeed nicely separated. Subsequent steps would then involve identification of marker genes, and more advanced downstream techniques, which are not part of the scope of this vignette. To read more on what can be done, the vignettes of the simpleSingleCell workflow package are an excellent place to start.

4.2 Exploring the dataset with iSEE

Once the processing steps above are done, we can call iSEE with the subsampled SingleCellExperiment object.

if (require(iSEE)) {
  iSEE(sce)
}

4.3 Saving the processed object

You can save the sce object to a serialized R object with

destination <- "where/to/store/the/processed/data.rds"
saveRDS(sce, file = destination)

The object can be read into a new R session with readRDS(destination), provided the HDF5 file remains in its original location (conveniently stored in the default location of ExperimentHub).