# (PART) Case studies {-}

# Human PBMCs (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

This performs an analysis of the public PBMC ID dataset generated by 10X Genomics [@zheng2017massively],
starting from the filtered count matrix.

## Data loading


```r
library(TENxPBMCData)
all.sce <- list(
    pbmc3k=TENxPBMCData('pbmc3k'),
    pbmc4k=TENxPBMCData('pbmc4k'),
    pbmc8k=TENxPBMCData('pbmc8k')
)
```

## Quality control


```r
unfiltered <- all.sce
```

Cell calling implicitly serves as a QC step to remove libraries with low total counts and number of detected genes.
Thus, we will only filter on the mitochondrial proportion.


```r
library(scater)
stats <- high.mito <- list()
for (n in names(all.sce)) {
    current <- all.sce[[n]]
    is.mito <- grep("MT", rowData(current)$Symbol_TENx)
    stats[[n]] <- perCellQCMetrics(current, subsets=list(Mito=is.mito))
    high.mito[[n]] <- isOutlier(stats[[n]]$subsets_Mito_percent, type="higher")
    all.sce[[n]] <- current[,!high.mito[[n]]]
}
```


```r
qcplots <- list()
for (n in names(all.sce)) {
    current <- unfiltered[[n]]
    colData(current) <- cbind(colData(current), stats[[n]])
    current$discard <- high.mito[[n]]
    qcplots[[n]] <- plotColData(current, x="sum", y="subsets_Mito_percent",
        colour_by="discard") + scale_x_log10()
}
do.call(gridExtra::grid.arrange, c(qcplots, ncol=3))
```

<div class="figure">
<img src="tenx-filtered-pbmc3k-4k-8k_files/figure-html/unref-pbmc-filtered-var-1.png" alt="Percentage of mitochondrial reads in each cell in each of the 10X PBMC datasets, compared to the total count. Each point represents a cell and is colored according to whether that cell was discarded." width="672" />
<p class="caption">(\#fig:unref-pbmc-filtered-var)Percentage of mitochondrial reads in each cell in each of the 10X PBMC datasets, compared to the total count. Each point represents a cell and is colored according to whether that cell was discarded.</p>
</div>


```r
lapply(high.mito, summary)
```

```
## $pbmc3k
##    Mode   FALSE    TRUE 
## logical    2609      91 
## 
## $pbmc4k
##    Mode   FALSE    TRUE 
## logical    4182     158 
## 
## $pbmc8k
##    Mode   FALSE    TRUE 
## logical    8157     224
```

## Normalization

We perform library size normalization, simply for convenience when dealing with file-backed matrices.


```r
all.sce <- lapply(all.sce, logNormCounts)
```


```r
lapply(all.sce, function(x) summary(sizeFactors(x)))
```

```
## $pbmc3k
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.234   0.748   0.926   1.000   1.157   6.604 
## 
## $pbmc4k
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.315   0.711   0.890   1.000   1.127  11.027 
## 
## $pbmc8k
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.296   0.704   0.877   1.000   1.118   6.794
```

## Variance modelling


```r
library(scran)
all.dec <- lapply(all.sce, modelGeneVar)
all.hvgs <- lapply(all.dec, getTopHVGs, prop=0.1)
```


```r
par(mfrow=c(1,3))
for (n in names(all.dec)) {
    curdec <- all.dec[[n]]
    plot(curdec$mean, curdec$total, pch=16, cex=0.5, main=n,
        xlab="Mean of log-expression", ylab="Variance of log-expression")
    curfit <- metadata(curdec)
    curve(curfit$trend(x), col='dodgerblue', add=TRUE, lwd=2)
}
```

<div class="figure">
<img src="tenx-filtered-pbmc3k-4k-8k_files/figure-html/unref-filtered-pbmc-variance-1.png" alt="Per-gene variance as a function of the mean for the log-expression values in each PBMC dataset. Each point represents a gene (black) with the mean-variance trend (blue) fitted to the variances." width="672" />
<p class="caption">(\#fig:unref-filtered-pbmc-variance)Per-gene variance as a function of the mean for the log-expression values in each PBMC dataset. Each point represents a gene (black) with the mean-variance trend (blue) fitted to the variances.</p>
</div>

## Dimensionality reduction

For various reasons, we will first analyze each PBMC dataset separately rather than merging them together.
We use randomized SVD, which is more efficient for file-backed matrices.


```r
library(BiocSingular)
set.seed(10000)
all.sce <- mapply(FUN=runPCA, x=all.sce, subset_row=all.hvgs, 
    MoreArgs=list(ncomponents=25, BSPARAM=RandomParam()), 
    SIMPLIFY=FALSE)

set.seed(100000)
all.sce <- lapply(all.sce, runTSNE, dimred="PCA")

set.seed(1000000)
all.sce <- lapply(all.sce, runUMAP, dimred="PCA")
```

## Clustering


```r
for (n in names(all.sce)) {
    g <- buildSNNGraph(all.sce[[n]], k=10, use.dimred='PCA')
    clust <- igraph::cluster_walktrap(g)$membership
    colLabels(all.sce[[n]])  <- factor(clust)
}
```


```r
lapply(all.sce, function(x) table(colLabels(x)))
```

```
## $pbmc3k
## 
##   1   2   3   4   5   6   7   8   9  10 
## 487 154 603 514  31 150 179 333 147  11 
## 
## $pbmc4k
## 
##    1    2    3    4    5    6    7    8    9   10   11   12   13 
##  497  185  569  786  373  232   44 1023   77  218   88   54   36 
## 
## $pbmc8k
## 
##    1    2    3    4    5    6    7    8    9   10   11   12   13   14   15   16 
## 1004  759 1073 1543  367  150  201 2067   59  154  244   67   76  285   20   15 
##   17   18 
##   64    9
```


```r
all.tsne <- list()
for (n in names(all.sce)) {
    all.tsne[[n]] <- plotTSNE(all.sce[[n]], colour_by="label") + ggtitle(n)
}
do.call(gridExtra::grid.arrange, c(all.tsne, list(ncol=2)))
```

<div class="figure">
<img src="tenx-filtered-pbmc3k-4k-8k_files/figure-html/unref-filtered-pbmc-tsne-1.png" alt="Obligatory $t$-SNE plots of each PBMC dataset, where each point represents a cell in the corresponding dataset and is colored according to the assigned cluster." width="672" />
<p class="caption">(\#fig:unref-filtered-pbmc-tsne)Obligatory $t$-SNE plots of each PBMC dataset, where each point represents a cell in the corresponding dataset and is colored according to the assigned cluster.</p>
</div>

## Data integration

With the per-dataset analyses out of the way, we will now repeat the analysis after merging together the three batches.


```r
# Intersecting the common genes.
universe <- Reduce(intersect, lapply(all.sce, rownames))
all.sce2 <- lapply(all.sce, "[", i=universe,)
all.dec2 <- lapply(all.dec, "[", i=universe,)

# Renormalizing to adjust for differences in depth.
library(batchelor)
normed.sce <- do.call(multiBatchNorm, all.sce2)

# Identifying a set of HVGs using stats from all batches.
combined.dec <- do.call(combineVar, all.dec2)
combined.hvg <- getTopHVGs(combined.dec, n=5000)

set.seed(1000101)
merged.pbmc <- do.call(fastMNN, c(normed.sce, 
    list(subset.row=combined.hvg, BSPARAM=RandomParam())))
```

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


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

```
##         pbmc3k    pbmc4k   pbmc8k
## [1,] 7.003e-03 3.126e-03 0.000000
## [2,] 7.137e-05 5.125e-05 0.003003
```

We proceed to clustering:


```r
g <- buildSNNGraph(merged.pbmc, use.dimred="corrected")
colLabels(merged.pbmc) <- factor(igraph::cluster_louvain(g)$membership)
table(colLabels(merged.pbmc), merged.pbmc$batch)
```

```
##     
##      pbmc3k pbmc4k pbmc8k
##   1     508    395    810
##   2     327    588   1125
##   3     185    122    218
##   4     150    180    294
##   5     175    340    581
##   6     292    542   1014
##   7     347    638   1202
##   8     436    760   1558
##   9      17     27    111
##   10    114    388    832
##   11     33    109    179
##   12     11     54    160
##   13     11      3      9
##   14      3     36     64
```

And visualization:


```r
set.seed(10101010)
merged.pbmc <- runTSNE(merged.pbmc, dimred="corrected")
gridExtra::grid.arrange(
    plotTSNE(merged.pbmc, colour_by="label", text_by="label", text_colour="red"),
    plotTSNE(merged.pbmc, colour_by="batch")
)
```

<div class="figure">
<img src="tenx-filtered-pbmc3k-4k-8k_files/figure-html/unref-filtered-pbmc-merged-tsne-1.png" alt="Obligatory $t$-SNE plots for the merged PBMC datasets, where each point represents a cell and is colored by cluster (top) or batch (bottom)." width="672" />
<p class="caption">(\#fig:unref-filtered-pbmc-merged-tsne)Obligatory $t$-SNE plots for the merged PBMC datasets, where each point represents a cell and is colored by cluster (top) or batch (bottom).</p>
</div>

## Session Info {-}

<button class="rebook-collapse">View session info</button>
<div class="rebook-content">
```
R version 4.3.2 Patched (2023-11-13 r85521)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 22.04.3 LTS

Matrix products: default
BLAS:   /home/biocbuild/bbs-3.18-bioc/R/lib/libRblas.so 
LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.0

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_GB              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       

time zone: America/New_York
tzcode source: system (glibc)

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

other attached packages:
 [1] batchelor_1.18.1            BiocSingular_1.18.0        
 [3] scran_1.30.2                scater_1.30.1              
 [5] ggplot2_3.4.4               scuttle_1.12.0             
 [7] TENxPBMCData_1.20.0         HDF5Array_1.30.1           
 [9] rhdf5_2.46.1                DelayedArray_0.28.0        
[11] SparseArray_1.2.4           S4Arrays_1.2.0             
[13] abind_1.4-5                 Matrix_1.6-5               
[15] SingleCellExperiment_1.24.0 SummarizedExperiment_1.32.0
[17] Biobase_2.62.0              GenomicRanges_1.54.1       
[19] GenomeInfoDb_1.38.6         IRanges_2.36.0             
[21] S4Vectors_0.40.2            BiocGenerics_0.48.1        
[23] MatrixGenerics_1.14.0       matrixStats_1.2.0          
[25] BiocStyle_2.30.0            rebook_1.12.0              

loaded via a namespace (and not attached):
  [1] jsonlite_1.8.8                CodeDepends_0.6.5            
  [3] magrittr_2.0.3                ggbeeswarm_0.7.2             
  [5] farver_2.1.1                  rmarkdown_2.25               
  [7] zlibbioc_1.48.0               vctrs_0.6.5                  
  [9] memoise_2.0.1                 DelayedMatrixStats_1.24.0    
 [11] RCurl_1.98-1.14               htmltools_0.5.7              
 [13] AnnotationHub_3.10.0          curl_5.2.0                   
 [15] BiocNeighbors_1.20.2          Rhdf5lib_1.24.2              
 [17] sass_0.4.8                    bslib_0.6.1                  
 [19] cachem_1.0.8                  ResidualMatrix_1.12.0        
 [21] igraph_2.0.1.1                mime_0.12                    
 [23] lifecycle_1.0.4               pkgconfig_2.0.3              
 [25] rsvd_1.0.5                    R6_2.5.1                     
 [27] fastmap_1.1.1                 GenomeInfoDbData_1.2.11      
 [29] shiny_1.8.0                   digest_0.6.34                
 [31] colorspace_2.1-0              AnnotationDbi_1.64.1         
 [33] dqrng_0.3.2                   irlba_2.3.5.1                
 [35] ExperimentHub_2.10.0          RSQLite_2.3.5                
 [37] beachmat_2.18.1               filelock_1.0.3               
 [39] labeling_0.4.3                fansi_1.0.6                  
 [41] httr_1.4.7                    compiler_4.3.2               
 [43] bit64_4.0.5                   withr_3.0.0                  
 [45] BiocParallel_1.36.0           viridis_0.6.5                
 [47] DBI_1.2.1                     highr_0.10                   
 [49] rappdirs_0.3.3                bluster_1.12.0               
 [51] tools_4.3.2                   vipor_0.4.7                  
 [53] beeswarm_0.4.0                interactiveDisplayBase_1.40.0
 [55] httpuv_1.6.14                 glue_1.7.0                   
 [57] rhdf5filters_1.14.1           promises_1.2.1               
 [59] grid_4.3.2                    Rtsne_0.17                   
 [61] cluster_2.1.6                 generics_0.1.3               
 [63] gtable_0.3.4                  metapod_1.10.1               
 [65] ScaledMatrix_1.10.0           utf8_1.2.4                   
 [67] XVector_0.42.0                RcppAnnoy_0.0.22             
 [69] ggrepel_0.9.5                 BiocVersion_3.18.1           
 [71] pillar_1.9.0                  limma_3.58.1                 
 [73] later_1.3.2                   dplyr_1.1.4                  
 [75] BiocFileCache_2.10.1          lattice_0.22-5               
 [77] FNN_1.1.4                     bit_4.0.5                    
 [79] tidyselect_1.2.0              locfit_1.5-9.8               
 [81] Biostrings_2.70.2             knitr_1.45                   
 [83] gridExtra_2.3                 bookdown_0.37                
 [85] edgeR_4.0.15                  xfun_0.42                    
 [87] statmod_1.5.0                 yaml_2.3.8                   
 [89] evaluate_0.23                 codetools_0.2-19             
 [91] tibble_3.2.1                  BiocManager_1.30.22          
 [93] graph_1.80.0                  cli_3.6.2                    
 [95] uwot_0.1.16                   xtable_1.8-4                 
 [97] munsell_0.5.0                 jquerylib_0.1.4              
 [99] Rcpp_1.0.12                   dir.expiry_1.10.0            
[101] dbplyr_2.4.0                  png_0.1-8                    
[103] XML_3.99-0.16.1               parallel_4.3.2               
[105] ellipsis_0.3.2                blob_1.2.4                   
[107] sparseMatrixStats_1.14.0      bitops_1.0-7                 
[109] viridisLite_0.4.2             scales_1.3.0                 
[111] purrr_1.0.2                   crayon_1.5.2                 
[113] rlang_1.1.3                   cowplot_1.1.3                
[115] KEGGREST_1.42.0              
```
</div>