# Unfiltered 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

Here, we describe a brief analysis of the peripheral blood mononuclear cell (PBMC) dataset from 10X Genomics [@zheng2017massively].
The data are publicly available from the [10X Genomics website](https://support.10xgenomics.com/single-cell-gene-expression/datasets/2.1.0/pbmc4k),
from which we download the raw gene/barcode count matrices, i.e., before cell calling from the _CellRanger_ pipeline.

## Data loading


```r
library(DropletTestFiles)
raw.path <- getTestFile("tenx-2.1.0-pbmc4k/1.0.0/raw.tar.gz")
out.path <- file.path(tempdir(), "pbmc4k")
untar(raw.path, exdir=out.path)

library(DropletUtils)
fname <- file.path(out.path, "raw_gene_bc_matrices/GRCh38")
sce.pbmc <- read10xCounts(fname, col.names=TRUE)
```


```r
library(scater)
rownames(sce.pbmc) <- uniquifyFeatureNames(
    rowData(sce.pbmc)$ID, rowData(sce.pbmc)$Symbol)

library(EnsDb.Hsapiens.v86)
location <- mapIds(EnsDb.Hsapiens.v86, keys=rowData(sce.pbmc)$ID, 
    column="SEQNAME", keytype="GENEID")
```

## Quality control

We perform cell detection using the `emptyDrops()` algorithm, as discussed in [Advanced Section 7.2](http://bioconductor.org/books/3.17/OSCA.advanced/droplet-processing.html#qc-droplets).


```r
set.seed(100)
e.out <- emptyDrops(counts(sce.pbmc))
sce.pbmc <- sce.pbmc[,which(e.out$FDR <= 0.001)]
```


```r
unfiltered <- sce.pbmc
```

We use a relaxed QC strategy and only remove cells with large mitochondrial proportions, using it as a proxy for cell damage.
This reduces the risk of removing cell types with low RNA content, especially in a heterogeneous PBMC population with many different cell types.


```r
stats <- perCellQCMetrics(sce.pbmc, subsets=list(Mito=which(location=="MT")))
high.mito <- isOutlier(stats$subsets_Mito_percent, type="higher")
sce.pbmc <- sce.pbmc[,!high.mito]
```


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

```
##    Mode   FALSE    TRUE 
## logical    3985     315
```


```r
colData(unfiltered) <- cbind(colData(unfiltered), stats)
unfiltered$discard <- high.mito

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

<div class="figure">
<img src="tenx-unfiltered-pbmc4k_files/figure-html/unref-unfiltered-pbmc-qc-1.png" alt="Distribution of various QC metrics in the PBMC dataset after cell calling. Each point is a cell and is colored according to whether it was discarded by the mitochondrial filter." width="672" />
<p class="caption">(\#fig:unref-unfiltered-pbmc-qc)Distribution of various QC metrics in the PBMC dataset after cell calling. Each point is a cell and is colored according to whether it was discarded by the mitochondrial filter.</p>
</div>


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

<div class="figure">
<img src="tenx-unfiltered-pbmc4k_files/figure-html/unref-unfiltered-pbmc-mito-1.png" alt="Proportion of mitochondrial reads in each cell of the PBMC dataset compared to its total count." width="672" />
<p class="caption">(\#fig:unref-unfiltered-pbmc-mito)Proportion of mitochondrial reads in each cell of the PBMC dataset compared to its total count.</p>
</div>

## Normalization


```r
library(scran)
set.seed(1000)
clusters <- quickCluster(sce.pbmc)
sce.pbmc <- computeSumFactors(sce.pbmc, cluster=clusters)
sce.pbmc <- logNormCounts(sce.pbmc)
```


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

```
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.007   0.712   0.875   1.000   1.099  12.254
```


```r
plot(librarySizeFactors(sce.pbmc), sizeFactors(sce.pbmc), pch=16,
    xlab="Library size factors", ylab="Deconvolution factors", log="xy")
```

<div class="figure">
<img src="tenx-unfiltered-pbmc4k_files/figure-html/unref-unfiltered-pbmc-norm-1.png" alt="Relationship between the library size factors and the deconvolution size factors in the PBMC dataset." width="672" />
<p class="caption">(\#fig:unref-unfiltered-pbmc-norm)Relationship between the library size factors and the deconvolution size factors in the PBMC dataset.</p>
</div>

## Variance modelling


```r
set.seed(1001)
dec.pbmc <- modelGeneVarByPoisson(sce.pbmc)
top.pbmc <- getTopHVGs(dec.pbmc, prop=0.1)
```


```r
plot(dec.pbmc$mean, dec.pbmc$total, pch=16, cex=0.5,
    xlab="Mean of log-expression", ylab="Variance of log-expression")
curfit <- metadata(dec.pbmc)
curve(curfit$trend(x), col='dodgerblue', add=TRUE, lwd=2)
```

<div class="figure">
<img src="tenx-unfiltered-pbmc4k_files/figure-html/unref-unfiltered-pbmc-var-1.png" alt="Per-gene variance as a function of the mean for the log-expression values in the PBMC dataset. Each point represents a gene (black) with the mean-variance trend (blue) fitted to simulated Poisson counts." width="672" />
<p class="caption">(\#fig:unref-unfiltered-pbmc-var)Per-gene variance as a function of the mean for the log-expression values in the PBMC dataset. Each point represents a gene (black) with the mean-variance trend (blue) fitted to simulated Poisson counts.</p>
</div>

## Dimensionality reduction


```r
set.seed(10000)
sce.pbmc <- denoisePCA(sce.pbmc, subset.row=top.pbmc, technical=dec.pbmc)

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

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

We verify that a reasonable number of PCs is retained.


```r
ncol(reducedDim(sce.pbmc, "PCA"))
```

```
## [1] 9
```

## Clustering


```r
g <- buildSNNGraph(sce.pbmc, k=10, use.dimred = 'PCA')
clust <- igraph::cluster_walktrap(g)$membership
colLabels(sce.pbmc) <- factor(clust)
```


```r
table(colLabels(sce.pbmc))
```

```
## 
##   1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16 
## 205 508 541  56 374 125  46 432 302 867  47 155 166  61  84  16
```


```r
plotTSNE(sce.pbmc, colour_by="label")
```

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

## Interpretation




```r
markers <- findMarkers(sce.pbmc, pval.type="some", direction="up")
```

We examine the markers for cluster 8 in more detail.
High expression of _CD14_, _CD68_ and _MNDA_ combined with low expression of _CD16_ suggests that this cluster contains monocytes,
compared to macrophages in cluster 15 (Figure \@ref(fig:unref-mono-pbmc-markers)).


```r
marker.set <- markers[["8"]]
as.data.frame(marker.set[1:30,1:3])
```

```
##                  p.value        FDR summary.logFC
## CSTA          7.171e-222 2.016e-217        2.4179
## MNDA          1.197e-221 2.016e-217        2.6615
## FCN1          2.376e-213 2.669e-209        2.6381
## S100A12       4.393e-212 3.701e-208        3.0809
## VCAN          1.711e-199 1.153e-195        2.2604
## TYMP          1.174e-154 6.590e-151        2.0238
## AIF1          3.674e-149 1.768e-145        2.4604
## LGALS2        4.005e-137 1.687e-133        1.8928
## MS4A6A        5.640e-134 2.111e-130        1.5457
## FGL2          2.045e-124 6.889e-121        1.3859
## RP11-1143G9.4 6.892e-122 2.111e-118        2.8042
## AP1S2         1.786e-112 5.015e-109        1.7704
## CD14          1.195e-110 3.098e-107        1.4260
## CFD           6.870e-109 1.654e-105        1.3560
## GPX1          9.049e-107 2.033e-103        2.4014
## TNFSF13B       3.920e-95  8.256e-92        1.1151
## KLF4           3.310e-94  6.560e-91        1.2049
## GRN            4.801e-91  8.987e-88        1.3815
## NAMPT          2.490e-90  4.415e-87        1.1439
## CLEC7A         7.736e-88  1.303e-84        1.0616
## S100A8         3.125e-84  5.014e-81        4.8052
## SERPINA1       1.580e-82  2.420e-79        1.3843
## CD36           8.018e-79  1.175e-75        1.0538
## MPEG1          8.482e-79  1.191e-75        0.9778
## CD68           5.119e-78  6.899e-75        0.9481
## CYBB           1.201e-77  1.556e-74        1.0300
## S100A11        1.175e-72  1.466e-69        1.8962
## RBP7           2.467e-71  2.969e-68        0.9666
## BLVRB          3.763e-71  4.372e-68        0.9701
## CD302          9.859e-71  1.107e-67        0.8792
```




```r
plotExpression(sce.pbmc, features=c("CD14", "CD68",
    "MNDA", "FCGR3A"), x="label", colour_by="label")
```

<div class="figure">
<img src="tenx-unfiltered-pbmc4k_files/figure-html/unref-mono-pbmc-markers-1.png" alt="Distribution of expression values for monocyte and macrophage markers across clusters in the PBMC dataset." width="672" />
<p class="caption">(\#fig:unref-mono-pbmc-markers)Distribution of expression values for monocyte and macrophage markers across clusters in the PBMC dataset.</p>
</div>

## Session Info {-}

<button class="rebook-collapse">View session info</button>
<div class="rebook-content">
```
R version 4.3.0 RC (2023-04-13 r84269)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 22.04.2 LTS

Matrix products: default
BLAS:   /home/biocbuild/bbs-3.17-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] scran_1.28.0                EnsDb.Hsapiens.v86_2.99.0  
 [3] ensembldb_2.24.0            AnnotationFilter_1.24.0    
 [5] GenomicFeatures_1.52.0      AnnotationDbi_1.62.0       
 [7] scater_1.28.0               ggplot2_3.4.2              
 [9] scuttle_1.10.0              DropletUtils_1.20.0        
[11] SingleCellExperiment_1.22.0 SummarizedExperiment_1.30.0
[13] Biobase_2.60.0              GenomicRanges_1.52.0       
[15] GenomeInfoDb_1.36.0         IRanges_2.34.0             
[17] S4Vectors_0.38.0            BiocGenerics_0.46.0        
[19] MatrixGenerics_1.12.0       matrixStats_0.63.0         
[21] DropletTestFiles_1.9.0      BiocStyle_2.28.0           
[23] rebook_1.10.0              

loaded via a namespace (and not attached):
  [1] later_1.3.0                   BiocIO_1.10.0                
  [3] bitops_1.0-7                  filelock_1.0.2               
  [5] tibble_3.2.1                  R.oo_1.25.0                  
  [7] CodeDepends_0.6.5             graph_1.78.0                 
  [9] XML_3.99-0.14                 lifecycle_1.0.3              
 [11] edgeR_3.42.0                  lattice_0.21-8               
 [13] magrittr_2.0.3                limma_3.56.0                 
 [15] sass_0.4.5                    rmarkdown_2.21               
 [17] jquerylib_0.1.4               yaml_2.3.7                   
 [19] metapod_1.8.0                 httpuv_1.6.9                 
 [21] cowplot_1.1.1                 DBI_1.1.3                    
 [23] zlibbioc_1.46.0               Rtsne_0.16                   
 [25] purrr_1.0.1                   R.utils_2.12.2               
 [27] RCurl_1.98-1.12               rappdirs_0.3.3               
 [29] GenomeInfoDbData_1.2.10       ggrepel_0.9.3                
 [31] irlba_2.3.5.1                 dqrng_0.3.0                  
 [33] DelayedMatrixStats_1.22.0     codetools_0.2-19             
 [35] DelayedArray_0.26.0           xml2_1.3.3                   
 [37] tidyselect_1.2.0              farver_2.1.1                 
 [39] ScaledMatrix_1.8.0            viridis_0.6.2                
 [41] BiocFileCache_2.8.0           GenomicAlignments_1.36.0     
 [43] jsonlite_1.8.4                BiocNeighbors_1.18.0         
 [45] ellipsis_0.3.2                tools_4.3.0                  
 [47] progress_1.2.2                Rcpp_1.0.10                  
 [49] glue_1.6.2                    gridExtra_2.3                
 [51] xfun_0.39                     dplyr_1.1.2                  
 [53] HDF5Array_1.28.0              withr_2.5.0                  
 [55] BiocManager_1.30.20           fastmap_1.1.1                
 [57] rhdf5filters_1.12.0           bluster_1.10.0               
 [59] fansi_1.0.4                   digest_0.6.31                
 [61] rsvd_1.0.5                    R6_2.5.1                     
 [63] mime_0.12                     colorspace_2.1-0             
 [65] biomaRt_2.56.0                RSQLite_2.3.1                
 [67] R.methodsS3_1.8.2             utf8_1.2.3                   
 [69] generics_0.1.3                rtracklayer_1.60.0           
 [71] FNN_1.1.3.2                   prettyunits_1.1.1            
 [73] httr_1.4.5                    uwot_0.1.14                  
 [75] pkgconfig_2.0.3               gtable_0.3.3                 
 [77] blob_1.2.4                    XVector_0.40.0               
 [79] htmltools_0.5.5               bookdown_0.33                
 [81] ProtGenerics_1.32.0           scales_1.2.1                 
 [83] png_0.1-8                     knitr_1.42                   
 [85] rjson_0.2.21                  curl_5.0.0                   
 [87] cachem_1.0.7                  rhdf5_2.44.0                 
 [89] stringr_1.5.0                 BiocVersion_3.17.1           
 [91] parallel_4.3.0                vipor_0.4.5                  
 [93] restfulr_0.0.15               pillar_1.9.0                 
 [95] grid_4.3.0                    vctrs_0.6.2                  
 [97] promises_1.2.0.1              BiocSingular_1.16.0          
 [99] dbplyr_2.3.2                  beachmat_2.16.0              
[101] xtable_1.8-4                  cluster_2.1.4                
[103] beeswarm_0.4.0                evaluate_0.20                
[105] cli_3.6.1                     locfit_1.5-9.7               
[107] compiler_4.3.0                Rsamtools_2.16.0             
[109] rlang_1.1.0                   crayon_1.5.2                 
[111] labeling_0.4.2                ggbeeswarm_0.7.1             
[113] stringi_1.7.12                viridisLite_0.4.1            
[115] BiocParallel_1.34.0           munsell_0.5.0                
[117] Biostrings_2.68.0             lazyeval_0.2.2               
[119] Matrix_1.5-4                  dir.expiry_1.8.0             
[121] ExperimentHub_2.8.0           hms_1.1.3                    
[123] sparseMatrixStats_1.12.0      bit64_4.0.5                  
[125] Rhdf5lib_1.22.0               KEGGREST_1.40.0              
[127] statmod_1.5.0                 shiny_1.7.4                  
[129] interactiveDisplayBase_1.38.0 highr_0.10                   
[131] AnnotationHub_3.8.0           igraph_1.4.2                 
[133] memoise_2.0.1                 bslib_0.4.2                  
[135] bit_4.0.5                    
```
</div>