# 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.18/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.1 (2023-06-16)
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] scran_1.30.0                EnsDb.Hsapiens.v86_2.99.0  
 [3] ensembldb_2.26.0            AnnotationFilter_1.26.0    
 [5] GenomicFeatures_1.54.0      AnnotationDbi_1.64.0       
 [7] scater_1.30.0               ggplot2_3.4.4              
 [9] scuttle_1.12.0              Matrix_1.6-1.1             
[11] DropletUtils_1.22.0         SingleCellExperiment_1.24.0
[13] SummarizedExperiment_1.32.0 Biobase_2.62.0             
[15] GenomicRanges_1.54.0        GenomeInfoDb_1.38.0        
[17] IRanges_2.36.0              S4Vectors_0.40.0           
[19] BiocGenerics_0.48.0         MatrixGenerics_1.14.0      
[21] matrixStats_1.0.0           DropletTestFiles_1.11.0    
[23] BiocStyle_2.30.0            rebook_1.12.0              

loaded via a namespace (and not attached):
  [1] later_1.3.1                   BiocIO_1.12.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.80.0                 
  [9] XML_3.99-0.14                 lifecycle_1.0.3              
 [11] edgeR_4.0.0                   lattice_0.22-5               
 [13] magrittr_2.0.3                limma_3.58.0                 
 [15] sass_0.4.7                    rmarkdown_2.25               
 [17] jquerylib_0.1.4               yaml_2.3.7                   
 [19] metapod_1.10.0                httpuv_1.6.12                
 [21] cowplot_1.1.1                 DBI_1.1.3                    
 [23] abind_1.4-5                   zlibbioc_1.48.0              
 [25] Rtsne_0.16                    purrr_1.0.2                  
 [27] R.utils_2.12.2                RCurl_1.98-1.12              
 [29] rappdirs_0.3.3                GenomeInfoDbData_1.2.11      
 [31] ggrepel_0.9.4                 irlba_2.3.5.1                
 [33] dqrng_0.3.1                   DelayedMatrixStats_1.24.0    
 [35] codetools_0.2-19              DelayedArray_0.28.0          
 [37] xml2_1.3.5                    tidyselect_1.2.0             
 [39] farver_2.1.1                  ScaledMatrix_1.10.0          
 [41] viridis_0.6.4                 BiocFileCache_2.10.0         
 [43] GenomicAlignments_1.38.0      jsonlite_1.8.7               
 [45] BiocNeighbors_1.20.0          ellipsis_0.3.2               
 [47] tools_4.3.1                   progress_1.2.2               
 [49] Rcpp_1.0.11                   glue_1.6.2                   
 [51] gridExtra_2.3                 SparseArray_1.2.0            
 [53] xfun_0.40                     dplyr_1.1.3                  
 [55] HDF5Array_1.30.0              withr_2.5.1                  
 [57] BiocManager_1.30.22           fastmap_1.1.1                
 [59] rhdf5filters_1.14.0           bluster_1.12.0               
 [61] fansi_1.0.5                   digest_0.6.33                
 [63] rsvd_1.0.5                    R6_2.5.1                     
 [65] mime_0.12                     colorspace_2.1-0             
 [67] biomaRt_2.58.0                RSQLite_2.3.1                
 [69] R.methodsS3_1.8.2             utf8_1.2.4                   
 [71] generics_0.1.3                FNN_1.1.3.2                  
 [73] rtracklayer_1.62.0            prettyunits_1.2.0            
 [75] httr_1.4.7                    S4Arrays_1.2.0               
 [77] uwot_0.1.16                   pkgconfig_2.0.3              
 [79] gtable_0.3.4                  blob_1.2.4                   
 [81] XVector_0.42.0                htmltools_0.5.6.1            
 [83] bookdown_0.36                 ProtGenerics_1.34.0          
 [85] scales_1.2.1                  png_0.1-8                    
 [87] knitr_1.44                    rstudioapi_0.15.0            
 [89] rjson_0.2.21                  curl_5.1.0                   
 [91] cachem_1.0.8                  rhdf5_2.46.0                 
 [93] stringr_1.5.0                 BiocVersion_3.18.0           
 [95] parallel_4.3.1                vipor_0.4.5                  
 [97] restfulr_0.0.15               pillar_1.9.0                 
 [99] grid_4.3.1                    vctrs_0.6.4                  
[101] promises_1.2.1                BiocSingular_1.18.0          
[103] dbplyr_2.3.4                  beachmat_2.18.0              
[105] xtable_1.8-4                  cluster_2.1.4                
[107] beeswarm_0.4.0                evaluate_0.22                
[109] cli_3.6.1                     locfit_1.5-9.8               
[111] compiler_4.3.1                Rsamtools_2.18.0             
[113] rlang_1.1.1                   crayon_1.5.2                 
[115] labeling_0.4.3                ggbeeswarm_0.7.2             
[117] stringi_1.7.12                viridisLite_0.4.2            
[119] BiocParallel_1.36.0           munsell_0.5.0                
[121] Biostrings_2.70.0             lazyeval_0.2.2               
[123] dir.expiry_1.10.0             ExperimentHub_2.10.0         
[125] hms_1.1.3                     sparseMatrixStats_1.14.0     
[127] bit64_4.0.5                   Rhdf5lib_1.24.0              
[129] KEGGREST_1.42.0               statmod_1.5.0                
[131] shiny_1.7.5.1                 interactiveDisplayBase_1.40.0
[133] AnnotationHub_3.10.0          igraph_1.5.1                 
[135] memoise_2.0.1                 bslib_0.5.1                  
[137] bit_4.0.5                    
```
</div>