# HCA human bone marrow (10X Genomics) 

<script>
document.addEventListener("click", function (event) {
    if (event.target.classList.contains("rebook-collapse")) {
        event.target.classList.toggle("active");
        var content = event.target.nextElementSibling;
        if (content.style.display === "block") {
            content.style.display = "none";
        } else {
            content.style.display = "block";
        }
    }
})
</script>

<style>
.rebook-collapse {
  background-color: #eee;
  color: #444;
  cursor: pointer;
  padding: 18px;
  width: 100%;
  border: none;
  text-align: left;
  outline: none;
  font-size: 15px;
}

.rebook-content {
  padding: 0 18px;
  display: none;
  overflow: hidden;
  background-color: #f1f1f1;
}
</style>

## Introduction

Here, we use an example dataset from the [Human Cell Atlas immune cell profiling project on bone marrow](https://preview.data.humancellatlas.org), which contains scRNA-seq data for 380,000 cells generated using the 10X Genomics technology.
This is a fairly big dataset that represents a good use case for the techniques in [Advanced Chapter 14](http://bioconductor.org/books/3.13/OSCA.advanced/dealing-with-big-data.html#dealing-with-big-data).

## Data loading

This dataset is loaded via the *[HCAData](https://bioconductor.org/packages/3.13/HCAData)* package, which provides a ready-to-use `SingleCellExperiment` object. 


```r
library(HCAData)
sce.bone <- HCAData('ica_bone_marrow', as.sparse=TRUE)
sce.bone$Donor <- sub("_.*", "", sce.bone$Barcode)
```

We use symbols in place of IDs for easier interpretation later.


```r
library(EnsDb.Hsapiens.v86)
rowData(sce.bone)$Chr <- mapIds(EnsDb.Hsapiens.v86, keys=rownames(sce.bone),
    column="SEQNAME", keytype="GENEID")

library(scater)
rownames(sce.bone) <- uniquifyFeatureNames(rowData(sce.bone)$ID,
    names = rowData(sce.bone)$Symbol)
```

## Quality control

Cell calling was not performed (see [here](https://s3.amazonaws.com/preview-ica-expression-data/Brief+ICA+Read+Me.pdf)) so we will perform QC using all metrics and block on the donor of origin during outlier detection.
We perform the calculation across multiple cores to speed things up.


```r
library(BiocParallel)
bpp <- MulticoreParam(8)
sce.bone <- unfiltered <- addPerCellQC(sce.bone, BPPARAM=bpp,
    subsets=list(Mito=which(rowData(sce.bone)$Chr=="MT")))

qc <- quickPerCellQC(colData(sce.bone), batch=sce.bone$Donor,
    sub.fields="subsets_Mito_percent")
sce.bone <- sce.bone[,!qc$discard]
```


```r
unfiltered$discard <- qc$discard

gridExtra::grid.arrange(
    plotColData(unfiltered, x="Donor", y="sum", colour_by="discard") +
        scale_y_log10() + ggtitle("Total count"),
    plotColData(unfiltered, x="Donor", y="detected", colour_by="discard") +
        scale_y_log10() + ggtitle("Detected features"),
    plotColData(unfiltered, x="Donor", y="subsets_Mito_percent",
        colour_by="discard") + ggtitle("Mito percent"),
    ncol=2
)
```

<div class="figure">
<img src="hca-bone-marrow_files/figure-html/unref-hca-bone-qc-1.png" alt="Distribution of QC metrics in the HCA bone marrow dataset. Each point represents a cell and is colored according to whether it was discarded." width="672" />
<p class="caption">(\#fig:unref-hca-bone-qc)Distribution of QC metrics in the HCA bone marrow dataset. Each point represents a cell and is colored according to whether it was discarded.</p>
</div>


```r
plotColData(unfiltered, x="sum", y="subsets_Mito_percent", 
    colour_by="discard") + scale_x_log10()
```

<div class="figure">
<img src="hca-bone-marrow_files/figure-html/unref-hca-bone-mito-1.png" alt="Percentage of mitochondrial reads in each cell in the HCA bone marrow dataset compared to its total count. Each point represents a cell and is colored according to whether that cell was discarded." width="672" />
<p class="caption">(\#fig:unref-hca-bone-mito)Percentage of mitochondrial reads in each cell in the HCA bone marrow dataset compared to its total count. Each point represents a cell and is colored according to whether that cell was discarded.</p>
</div>

## Normalization

For a minor speed-up, we use already-computed library sizes rather than re-computing them from the column sums.


```r
sce.bone <- logNormCounts(sce.bone, size_factors = sce.bone$sum)
```


```r
summary(sizeFactors(sce.bone))
```

```
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.05    0.47    0.65    1.00    0.89   42.38
```

## Variance modeling

We block on the donor of origin to mitigate batch effects during HVG selection.
We select a larger number of HVGs to capture any batch-specific variation that might be present.


```r
library(scran)
set.seed(1010010101)
dec.bone <- modelGeneVarByPoisson(sce.bone, 
    block=sce.bone$Donor, BPPARAM=bpp)
top.bone <- getTopHVGs(dec.bone, n=5000)
```


```r
par(mfrow=c(4,2))
blocked.stats <- dec.bone$per.block
for (i in colnames(blocked.stats)) {
    current <- blocked.stats[[i]]
    plot(current$mean, current$total, main=i, pch=16, cex=0.5,
        xlab="Mean of log-expression", ylab="Variance of log-expression")
    curfit <- metadata(current)
    curve(curfit$trend(x), col='dodgerblue', add=TRUE, lwd=2)
}
```

<div class="figure">
<img src="hca-bone-marrow_files/figure-html/unref-hca-bone-var-1.png" alt="Per-gene variance as a function of the mean for the log-expression values in the HCA bone marrow dataset. Each point represents a gene (black) with the mean-variance trend (blue) fitted to the variances." width="576" />
<p class="caption">(\#fig:unref-hca-bone-var)Per-gene variance as a function of the mean for the log-expression values in the HCA bone marrow dataset. Each point represents a gene (black) with the mean-variance trend (blue) fitted to the variances.</p>
</div>

## Data integration

Here we use multiple cores, randomized SVD and approximate nearest-neighbor detection to speed up this step.


```r
library(batchelor)
library(BiocNeighbors)

set.seed(1010001)
merged.bone <- fastMNN(sce.bone, batch = sce.bone$Donor, subset.row = top.bone,
     BSPARAM=BiocSingular::RandomParam(deferred = TRUE), 
     BNPARAM=AnnoyParam(),
     BPPARAM=bpp)

reducedDim(sce.bone, 'MNN') <- reducedDim(merged.bone, 'corrected')
```

We use the percentage of variance lost as a diagnostic measure:


```r
metadata(merged.bone)$merge.info$lost.var
```

```
##      MantonBM1 MantonBM2 MantonBM3 MantonBM4 MantonBM5 MantonBM6 MantonBM7
## [1,]  0.006922  0.006392  0.000000  0.000000  0.000000  0.000000  0.000000
## [2,]  0.006380  0.006863  0.023049  0.000000  0.000000  0.000000  0.000000
## [3,]  0.005068  0.003084  0.005178  0.019496  0.000000  0.000000  0.000000
## [4,]  0.002009  0.001891  0.001901  0.001786  0.023105  0.000000  0.000000
## [5,]  0.002452  0.002003  0.001770  0.002926  0.002646  0.023852  0.000000
## [6,]  0.003167  0.003222  0.003169  0.002636  0.003362  0.003442  0.024650
## [7,]  0.001968  0.001701  0.002441  0.002045  0.001585  0.002312  0.002003
##      MantonBM8
## [1,]   0.00000
## [2,]   0.00000
## [3,]   0.00000
## [4,]   0.00000
## [5,]   0.00000
## [6,]   0.00000
## [7,]   0.03216
```

## Dimensionality reduction

We set `external_neighbors=TRUE` to replace the internal nearest neighbor search in the UMAP implementation with our parallelized approximate search.
We also set the number of threads to be used in the UMAP iterations.


```r
set.seed(01010100)
sce.bone <- runUMAP(sce.bone, dimred="MNN",
    external_neighbors=TRUE, 
    BNPARAM=AnnoyParam(),
    BPPARAM=bpp,
    n_threads=bpnworkers(bpp))
```

## Clustering

Graph-based clustering generates an excessively large intermediate graph so we will instead use a two-step approach with $k$-means.
We generate 1000 small clusters that are subsequently aggregated into more interpretable groups with a graph-based method.
If more resolution is required, we can increase `centers` in addition to using a lower `k` during graph construction.


```r
library(bluster)

set.seed(1000)
colLabels(sce.bone) <- clusterRows(reducedDim(sce.bone, "MNN"),
    TwoStepParam(KmeansParam(centers=1000), NNGraphParam(k=5)))

table(colLabels(sce.bone))
```

```
## 
##     1     2     3     4     5     6     7     8     9    10    11    12    13 
## 20331 11161 55464 47426 15731 10581 64721 26493 18703 15043 17097  4992  3157 
##    14    15 
##  3403  2422
```

We observe mostly balanced contributions from different samples to each cluster (Figure \@ref(fig:unref-hca-bone-ab)), consistent with the expectation that all samples are replicates from different donors.


```r
tab <- table(Cluster=colLabels(sce.bone), Donor=sce.bone$Donor)
library(pheatmap)
pheatmap(log10(tab+10), color=viridis::viridis(100))
```

<div class="figure">
<img src="hca-bone-marrow_files/figure-html/unref-hca-bone-ab-1.png" alt="Heatmap of log~10~-number of cells in each cluster (row) from each sample (column)." width="672" />
<p class="caption">(\#fig:unref-hca-bone-ab)Heatmap of log~10~-number of cells in each cluster (row) from each sample (column).</p>
</div>




```r
# TODO: add scrambling option in scater's plotting functions.
scrambled <- sample(ncol(sce.bone))

gridExtra::grid.arrange(
    plotUMAP(sce.bone, colour_by="label", text_by="label"),
    plotUMAP(sce.bone[,scrambled], colour_by="Donor")
)
```

<div class="figure">
<img src="hca-bone-marrow_files/figure-html/unref-hca-bone-umap-1.png" alt="UMAP plots of the HCA bone marrow dataset after merging. Each point represents a cell and is colored according to the assigned cluster (top) or the donor of origin (bottom)." width="672" />
<p class="caption">(\#fig:unref-hca-bone-umap)UMAP plots of the HCA bone marrow dataset after merging. Each point represents a cell and is colored according to the assigned cluster (top) or the donor of origin (bottom).</p>
</div>

## Differential expression

We identify marker genes for each cluster while blocking on the donor.


```r
markers.bone <- findMarkers(sce.bone, block = sce.bone$Donor, 
    direction = 'up', lfc = 1, BPPARAM=bpp)
```

We visualize the top markers for a randomly chosen cluster using a "dot plot" in Figure \@ref(fig:unref-hca-bone-dotplot).
The presence of upregulated genes like _LYZ_, _S100A8_ and _VCAN_ is consistent with a monocyte identity for this cluster.




```r
top.markers <- markers.bone[["4"]]
best <- top.markers[top.markers$Top <= 10,]
lfcs <- getMarkerEffects(best)

library(pheatmap)
pheatmap(lfcs, breaks=seq(-5, 5, length.out=101))
```

<div class="figure">
<img src="hca-bone-marrow_files/figure-html/unref-hca-bone-dotplot-1.png" alt="Heatmap of log~2~-fold changes for the top marker genes (rows) of cluster 4 compared to all other clusters (columns)." width="672" />
<p class="caption">(\#fig:unref-hca-bone-dotplot)Heatmap of log~2~-fold changes for the top marker genes (rows) of cluster 4 compared to all other clusters (columns).</p>
</div>



## Cell type classification

We perform automated cell type classification using a reference dataset to annotate each cluster based on its pseudo-bulk profile. 
This is faster than the per-cell approaches described in Chapter \@ref(cell-type-annotation) at the cost of the resolution required to detect heterogeneity inside a cluster.
Nonetheless, it is often sufficient for a quick assignment of cluster identity, and indeed, cluster 4 is also identified as consisting of monocytes from this analysis.


```r
se.aggregated <- sumCountsAcrossCells(sce.bone, id=colLabels(sce.bone), BPPARAM=bpp)

library(celldex)
hpc <- HumanPrimaryCellAtlasData()

library(SingleR)
anno.single <- SingleR(se.aggregated, ref = hpc, labels = hpc$label.main,
    assay.type.test="sum")
anno.single
```

```
## DataFrame with 15 rows and 5 columns
##                             scores     first.labels      tuning.scores
##                           <matrix>      <character>        <DataFrame>
## 1   0.366050:0.741975:0.637800:...              GMP  0.676512:0.339887
## 2   0.399229:0.715314:0.628516:... Pro-B_cell_CD34+  0.845406:0.776315
## 3   0.325812:0.654869:0.570039:...          NK_cell  0.612492:0.483279
## 4   0.296455:0.742466:0.529404:... Pre-B_cell_CD34-  0.564310:0.254385
## 5   0.345773:0.565378:0.479722:...          T_cells  0.616848:0.122525
## ...                            ...              ...                ...
## 11  0.326882:0.646229:0.561060:...          NK_cell 0.546934:0.5190984
## 12  0.380710:0.684123:0.784540:...              MEP 0.375942:0.3754427
## 13  0.368546:0.652935:0.580330:...           B_cell 0.527988:0.3268898
## 14  0.294361:0.706019:0.527282:... Pre-B_cell_CD34- 0.600000:0.2447875
## 15  0.339786:0.687074:0.569933:... Pre-B_cell_CD34- 0.231710:0.0998799
##               labels    pruned.labels
##          <character>      <character>
## 1                GMP              GMP
## 2   Pro-B_cell_CD34+ Pro-B_cell_CD34+
## 3            T_cells          T_cells
## 4           Monocyte         Monocyte
## 5            T_cells          T_cells
## ...              ...              ...
## 11           T_cells          T_cells
## 12        BM & Prog.       BM & Prog.
## 13            B_cell               NA
## 14          Monocyte         Monocyte
## 15               GMP              GMP
```



<!-- ```{r saveRDS, include=FALSE} -->
<!-- ## TEMP FIX -->
<!-- saveRDS(sce.bone, 'sce.bone.rds') -->
<!-- saveRDS(markers.bone, 'markers.bone.rds') -->

<!-- sce.bone <- readRDS('sce.bone.rds') -->
<!-- markers.bone <- readRDS('markers.bone.rds') -->
<!-- ``` -->

<!--
## Gene set Enrichment Analysis

To interpret differential expression (DE) results, gene set enrichment analysis (GSEA) can be performed to annotate DE results. The first step is to acquire a mapping of ENTREZID to gene symbol.


```r
## Get mappings of ENTREZID to Symbol
library(org.Hs.eg.db)
keys_entrez <- keys(org.Hs.eg.db, 'ENTREZID')
mapping_es <- AnnotationDbi::select(org.Hs.eg.db,
                                    keys = keys_entrez,
                                    columns = c('ENTREZID', 'SYMBOL'),
                                    keytype = 'ENTREZID')
mapping_es$ENTREZID <- as.integer(mapping_es$ENTREZID)
```

The next is to parse the gene sets (here we use a "tidy" representation of MSigDB gene sets) and remap them.


```r
## Get pathways of interest - convert to list with symbol
## devtools::install_github('stephenturner/msigdf') # "tidy" MSigDB genesets
library(msigdf)
library(dplyr)
mdb <- dplyr::inner_join(msigdf.human, mapping_es,
                         by = c('entrez' = 'ENTREZID')) %>%
    dplyr::filter(collection == 'c7') %>%
    dplyr::select(-collection, -entrez) %>%
    dplyr::group_nest(geneset)
pathways <- purrr::map(mdb$data, function(x) {
    as.character(x$SYMBOL)
})
names(pathways) <- mdb$geneset
```

Finally, the DE results from the `findMarkers()` function are used to test for differential gene set enrichment using the `fgsea` package. Here we compare cluster 3 against cluster 15, which are at opposite ends of a UMAP group.


```r
## Get stats based on markers search
## Compare cluster 3 (Dendritic cells) and 15 (HSCs)
stats <- markers.bone[[3]]$logFC.15
names(stats) <- rownames(markers.bone[[3]])

## Run fast gene set enrichment analysis
## Returns a data.table object
library(fgsea)
gse <- fgsea(pathways = pathways, 
             stats = stats,
             minSize = 15,
             maxSize = 500,
             nperm = 5000)
```

The enrichment scores of the top up and down pathways can then be plotted as a table:


```r
## Get top up and down pathways based on NES
gse.up <- gse[order(gse$NES, decreasing = TRUE)[1:2], ]$pathway
gse.dn <- gse[order(gse$NES)[1:2], ]$pathway
topPathways <- c(gse.up, rev(gse.dn))

## Plot the GSEA table
plotGseaTable(pathways[topPathways], stats,
              gse, 
              gseaParam = 0.5)
```

Or select pathways can also be plotted as a more traditional running score plot:


```r
## Traditional GSEA plot
plotEnrichment(pathways[[topPathways[1]]],
               stats)
```

### Trajectory analysis

Slingshot relies on providing a set of clusters on which to construct a trajectory.


```r
library(slingshot)

## Clusters to attempt trajectory on:
slc.clusters <- c(15, 22, 11, 7, 1, 3)

## Subsample dataset for speed
slc.dat <- sce.bone[, sce.bone$cluster %in% slc.clusters]
slc.subdat <- slc.dat[, sample(ncol(slc.dat), 20000)]

## Run slingshot
slc <- slingshot(slc.subdat,
                 clusterLabels = 'cluster',
                 reducedDim = 'UMAP')
```

The calculated principal curve can then be visualized on the UMAP space.


```r
## Slingshot trajectory plot
library(RColorBrewer)
colors <- colorRampPalette(brewer.pal(11, 'Spectral')[-6])(100)
plotcol <- colors[cut(slc$slingPseudotime_1, breaks = 100)]
plot(reducedDims(slc)$UMAP, col = plotcol, pch=16, asp = 1)
lines(SlingshotDataSet(slc), lwd = 2, col = 'black')
```
-->

## Session Info {-}

<button class="rebook-collapse">View session info</button>
<div class="rebook-content">
```
R version 4.1.0 (2021-05-18)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 20.04.2 LTS

Matrix products: default
BLAS:   /home/biocbuild/bbs-3.13-bioc/R/lib/libRblas.so
LAPACK: /home/biocbuild/bbs-3.13-bioc/R/lib/libRlapack.so

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=C              
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] parallel  stats4    stats     graphics  grDevices utils     datasets 
[8] methods   base     

other attached packages:
 [1] SingleR_1.6.0               celldex_1.2.0              
 [3] pheatmap_1.0.12             bluster_1.2.0              
 [5] BiocNeighbors_1.10.0        batchelor_1.8.0            
 [7] scran_1.20.0                BiocParallel_1.26.0        
 [9] scater_1.20.0               ggplot2_3.3.3              
[11] scuttle_1.2.0               EnsDb.Hsapiens.v86_2.99.0  
[13] ensembldb_2.16.0            AnnotationFilter_1.16.0    
[15] GenomicFeatures_1.44.0      AnnotationDbi_1.54.0       
[17] rhdf5_2.36.0                HCAData_1.8.0              
[19] SingleCellExperiment_1.14.0 SummarizedExperiment_1.22.0
[21] Biobase_2.52.0              GenomicRanges_1.44.0       
[23] GenomeInfoDb_1.28.0         IRanges_2.26.0             
[25] S4Vectors_0.30.0            BiocGenerics_0.38.0        
[27] MatrixGenerics_1.4.0        matrixStats_0.58.0         
[29] BiocStyle_2.20.0            rebook_1.2.0               

loaded via a namespace (and not attached):
  [1] AnnotationHub_3.0.0           BiocFileCache_2.0.0          
  [3] igraph_1.2.6                  lazyeval_0.2.2               
  [5] digest_0.6.27                 htmltools_0.5.1.1            
  [7] viridis_0.6.1                 fansi_0.4.2                  
  [9] magrittr_2.0.1                memoise_2.0.0                
 [11] ScaledMatrix_1.0.0            cluster_2.1.2                
 [13] limma_3.48.0                  Biostrings_2.60.0            
 [15] prettyunits_1.1.1             colorspace_2.0-1             
 [17] blob_1.2.1                    rappdirs_0.3.3               
 [19] xfun_0.23                     dplyr_1.0.6                  
 [21] crayon_1.4.1                  RCurl_1.98-1.3               
 [23] jsonlite_1.7.2                graph_1.70.0                 
 [25] glue_1.4.2                    gtable_0.3.0                 
 [27] zlibbioc_1.38.0               XVector_0.32.0               
 [29] DelayedArray_0.18.0           BiocSingular_1.8.0           
 [31] Rhdf5lib_1.14.0               HDF5Array_1.20.0             
 [33] scales_1.1.1                  edgeR_3.34.0                 
 [35] DBI_1.1.1                     Rcpp_1.0.6                   
 [37] viridisLite_0.4.0             xtable_1.8-4                 
 [39] progress_1.2.2                dqrng_0.3.0                  
 [41] bit_4.0.4                     rsvd_1.0.5                   
 [43] ResidualMatrix_1.2.0          metapod_1.0.0                
 [45] httr_1.4.2                    RColorBrewer_1.1-2           
 [47] dir.expiry_1.0.0              ellipsis_0.3.2               
 [49] pkgconfig_2.0.3               XML_3.99-0.6                 
 [51] farver_2.1.0                  uwot_0.1.10                  
 [53] CodeDepends_0.6.5             sass_0.4.0                   
 [55] dbplyr_2.1.1                  locfit_1.5-9.4               
 [57] utf8_1.2.1                    labeling_0.4.2               
 [59] tidyselect_1.1.1              rlang_0.4.11                 
 [61] later_1.2.0                   munsell_0.5.0                
 [63] BiocVersion_3.13.1            tools_4.1.0                  
 [65] cachem_1.0.5                  generics_0.1.0               
 [67] RSQLite_2.2.7                 ExperimentHub_2.0.0          
 [69] evaluate_0.14                 stringr_1.4.0                
 [71] fastmap_1.1.0                 yaml_2.2.1                   
 [73] knitr_1.33                    bit64_4.0.5                  
 [75] purrr_0.3.4                   KEGGREST_1.32.0              
 [77] sparseMatrixStats_1.4.0       mime_0.10                    
 [79] biomaRt_2.48.0                compiler_4.1.0               
 [81] beeswarm_0.3.1                filelock_1.0.2               
 [83] curl_4.3.1                    png_0.1-7                    
 [85] interactiveDisplayBase_1.30.0 statmod_1.4.36               
 [87] tibble_3.1.2                  bslib_0.2.5.1                
 [89] stringi_1.6.2                 highr_0.9                    
 [91] RSpectra_0.16-0               lattice_0.20-44              
 [93] ProtGenerics_1.24.0           Matrix_1.3-3                 
 [95] vctrs_0.3.8                   pillar_1.6.1                 
 [97] lifecycle_1.0.0               rhdf5filters_1.4.0           
 [99] BiocManager_1.30.15           jquerylib_0.1.4              
[101] cowplot_1.1.1                 bitops_1.0-7                 
[103] irlba_2.3.3                   httpuv_1.6.1                 
[105] rtracklayer_1.52.0            R6_2.5.0                     
[107] BiocIO_1.2.0                  bookdown_0.22                
[109] promises_1.2.0.1              gridExtra_2.3                
[111] vipor_0.4.5                   codetools_0.2-18             
[113] assertthat_0.2.1              rjson_0.2.20                 
[115] withr_2.4.2                   GenomicAlignments_1.28.0     
[117] Rsamtools_2.8.0               GenomeInfoDbData_1.2.6       
[119] hms_1.1.0                     grid_4.1.0                   
[121] beachmat_2.8.0                rmarkdown_2.8                
[123] DelayedMatrixStats_1.14.0     shiny_1.6.0                  
[125] ggbeeswarm_0.6.0              restfulr_0.0.13              
```
</div>