# 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.16/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.2.1 (2022-06-23)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 20.04.5 LTS

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

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       

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

other attached packages:
 [1] scran_1.26.0                EnsDb.Hsapiens.v86_2.99.0  
 [3] ensembldb_2.22.0            AnnotationFilter_1.22.0    
 [5] GenomicFeatures_1.50.0      AnnotationDbi_1.60.0       
 [7] scater_1.26.0               ggplot2_3.3.6              
 [9] scuttle_1.8.0               DropletUtils_1.18.0        
[11] SingleCellExperiment_1.20.0 SummarizedExperiment_1.28.0
[13] Biobase_2.58.0              GenomicRanges_1.50.0       
[15] GenomeInfoDb_1.34.0         IRanges_2.32.0             
[17] S4Vectors_0.36.0            BiocGenerics_0.44.0        
[19] MatrixGenerics_1.10.0       matrixStats_0.62.0         
[21] DropletTestFiles_1.7.0      BiocStyle_2.26.0           
[23] rebook_1.8.0               

loaded via a namespace (and not attached):
  [1] AnnotationHub_3.6.0           BiocFileCache_2.6.0          
  [3] igraph_1.3.5                  lazyeval_0.2.2               
  [5] BiocParallel_1.32.0           digest_0.6.30                
  [7] htmltools_0.5.3               viridis_0.6.2                
  [9] fansi_1.0.3                   magrittr_2.0.3               
 [11] memoise_2.0.1                 ScaledMatrix_1.6.0           
 [13] cluster_2.1.4                 limma_3.54.0                 
 [15] Biostrings_2.66.0             R.utils_2.12.1               
 [17] prettyunits_1.1.1             colorspace_2.0-3             
 [19] blob_1.2.3                    rappdirs_0.3.3               
 [21] ggrepel_0.9.1                 xfun_0.34                    
 [23] dplyr_1.0.10                  crayon_1.5.2                 
 [25] RCurl_1.98-1.9                jsonlite_1.8.3               
 [27] graph_1.76.0                  glue_1.6.2                   
 [29] gtable_0.3.1                  zlibbioc_1.44.0              
 [31] XVector_0.38.0                DelayedArray_0.24.0          
 [33] BiocSingular_1.14.0           Rhdf5lib_1.20.0              
 [35] HDF5Array_1.26.0              scales_1.2.1                 
 [37] DBI_1.1.3                     edgeR_3.40.0                 
 [39] Rcpp_1.0.9                    viridisLite_0.4.1            
 [41] xtable_1.8-4                  progress_1.2.2               
 [43] dqrng_0.3.0                   bit_4.0.4                    
 [45] rsvd_1.0.5                    metapod_1.6.0                
 [47] httr_1.4.4                    FNN_1.1.3.1                  
 [49] dir.expiry_1.6.0              ellipsis_0.3.2               
 [51] farver_2.1.1                  pkgconfig_2.0.3              
 [53] XML_3.99-0.12                 R.methodsS3_1.8.2            
 [55] uwot_0.1.14                   CodeDepends_0.6.5            
 [57] sass_0.4.2                    dbplyr_2.2.1                 
 [59] locfit_1.5-9.6                utf8_1.2.2                   
 [61] labeling_0.4.2                tidyselect_1.2.0             
 [63] rlang_1.0.6                   later_1.3.0                  
 [65] munsell_0.5.0                 BiocVersion_3.16.0           
 [67] tools_4.2.1                   cachem_1.0.6                 
 [69] cli_3.4.1                     generics_0.1.3               
 [71] RSQLite_2.2.18                ExperimentHub_2.6.0          
 [73] evaluate_0.17                 stringr_1.4.1                
 [75] fastmap_1.1.0                 yaml_2.3.6                   
 [77] knitr_1.40                    bit64_4.0.5                  
 [79] purrr_0.3.5                   KEGGREST_1.38.0              
 [81] sparseMatrixStats_1.10.0      mime_0.12                    
 [83] R.oo_1.25.0                   xml2_1.3.3                   
 [85] biomaRt_2.54.0                compiler_4.2.1               
 [87] beeswarm_0.4.0                filelock_1.0.2               
 [89] curl_4.3.3                    png_0.1-7                    
 [91] interactiveDisplayBase_1.36.0 statmod_1.4.37               
 [93] tibble_3.1.8                  bslib_0.4.0                  
 [95] stringi_1.7.8                 highr_0.9                    
 [97] bluster_1.8.0                 lattice_0.20-45              
 [99] ProtGenerics_1.30.0           Matrix_1.5-1                 
[101] vctrs_0.5.0                   pillar_1.8.1                 
[103] lifecycle_1.0.3               rhdf5filters_1.10.0          
[105] BiocManager_1.30.19           jquerylib_0.1.4              
[107] BiocNeighbors_1.16.0          cowplot_1.1.1                
[109] bitops_1.0-7                  irlba_2.3.5.1                
[111] rtracklayer_1.58.0            httpuv_1.6.6                 
[113] R6_2.5.1                      BiocIO_1.8.0                 
[115] bookdown_0.29                 promises_1.2.0.1             
[117] gridExtra_2.3                 vipor_0.4.5                  
[119] codetools_0.2-18              assertthat_0.2.1             
[121] rhdf5_2.42.0                  rjson_0.2.21                 
[123] withr_2.5.0                   GenomicAlignments_1.34.0     
[125] Rsamtools_2.14.0              GenomeInfoDbData_1.2.9       
[127] parallel_4.2.1                hms_1.1.2                    
[129] grid_4.2.1                    beachmat_2.14.0              
[131] rmarkdown_2.17                DelayedMatrixStats_1.20.0    
[133] Rtsne_0.16                    shiny_1.7.3                  
[135] ggbeeswarm_0.6.0              restfulr_0.0.15              
```
</div>