# 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.13/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.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] scran_1.20.0                EnsDb.Hsapiens.v86_2.99.0  
 [3] ensembldb_2.16.0            AnnotationFilter_1.16.0    
 [5] GenomicFeatures_1.44.0      AnnotationDbi_1.54.0       
 [7] scater_1.20.0               ggplot2_3.3.3              
 [9] scuttle_1.2.0               DropletUtils_1.12.0        
[11] SingleCellExperiment_1.14.0 SummarizedExperiment_1.22.0
[13] Biobase_2.52.0              GenomicRanges_1.44.0       
[15] GenomeInfoDb_1.28.0         IRanges_2.26.0             
[17] S4Vectors_0.30.0            BiocGenerics_0.38.0        
[19] MatrixGenerics_1.4.0        matrixStats_0.58.0         
[21] DropletTestFiles_1.2.0      BiocStyle_2.20.0           
[23] 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] BiocParallel_1.26.0           digest_0.6.27                
  [7] htmltools_0.5.1.1             viridis_0.6.1                
  [9] fansi_0.4.2                   magrittr_2.0.1               
 [11] memoise_2.0.0                 ScaledMatrix_1.0.0           
 [13] cluster_2.1.2                 limma_3.48.0                 
 [15] Biostrings_2.60.0             R.utils_2.10.1               
 [17] prettyunits_1.1.1             colorspace_2.0-1             
 [19] blob_1.2.1                    rappdirs_0.3.3               
 [21] xfun_0.23                     dplyr_1.0.6                  
 [23] crayon_1.4.1                  RCurl_1.98-1.3               
 [25] jsonlite_1.7.2                graph_1.70.0                 
 [27] glue_1.4.2                    gtable_0.3.0                 
 [29] zlibbioc_1.38.0               XVector_0.32.0               
 [31] DelayedArray_0.18.0           BiocSingular_1.8.0           
 [33] Rhdf5lib_1.14.0               HDF5Array_1.20.0             
 [35] scales_1.1.1                  DBI_1.1.1                    
 [37] edgeR_3.34.0                  Rcpp_1.0.6                   
 [39] viridisLite_0.4.0             xtable_1.8-4                 
 [41] progress_1.2.2                dqrng_0.3.0                  
 [43] bit_4.0.4                     rsvd_1.0.5                   
 [45] metapod_1.0.0                 httr_1.4.2                   
 [47] FNN_1.1.3                     dir.expiry_1.0.0             
 [49] ellipsis_0.3.2                farver_2.1.0                 
 [51] pkgconfig_2.0.3               XML_3.99-0.6                 
 [53] R.methodsS3_1.8.1             uwot_0.1.10                  
 [55] CodeDepends_0.6.5             sass_0.4.0                   
 [57] dbplyr_2.1.1                  locfit_1.5-9.4               
 [59] utf8_1.2.1                    labeling_0.4.2               
 [61] tidyselect_1.1.1              rlang_0.4.11                 
 [63] later_1.2.0                   munsell_0.5.0                
 [65] BiocVersion_3.13.1            tools_4.1.0                  
 [67] cachem_1.0.5                  generics_0.1.0               
 [69] RSQLite_2.2.7                 ExperimentHub_2.0.0          
 [71] evaluate_0.14                 stringr_1.4.0                
 [73] fastmap_1.1.0                 yaml_2.2.1                   
 [75] knitr_1.33                    bit64_4.0.5                  
 [77] purrr_0.3.4                   KEGGREST_1.32.0              
 [79] sparseMatrixStats_1.4.0       mime_0.10                    
 [81] R.oo_1.24.0                   biomaRt_2.48.0               
 [83] compiler_4.1.0                beeswarm_0.3.1               
 [85] filelock_1.0.2                curl_4.3.1                   
 [87] png_0.1-7                     interactiveDisplayBase_1.30.0
 [89] statmod_1.4.36                tibble_3.1.2                 
 [91] bslib_0.2.5.1                 stringi_1.6.2                
 [93] highr_0.9                     RSpectra_0.16-0              
 [95] bluster_1.2.0                 lattice_0.20-44              
 [97] ProtGenerics_1.24.0           Matrix_1.3-3                 
 [99] vctrs_0.3.8                   pillar_1.6.1                 
[101] lifecycle_1.0.0               rhdf5filters_1.4.0           
[103] BiocManager_1.30.15           jquerylib_0.1.4              
[105] BiocNeighbors_1.10.0          cowplot_1.1.1                
[107] bitops_1.0-7                  irlba_2.3.3                  
[109] httpuv_1.6.1                  rtracklayer_1.52.0           
[111] R6_2.5.0                      BiocIO_1.2.0                 
[113] bookdown_0.22                 promises_1.2.0.1             
[115] gridExtra_2.3                 vipor_0.4.5                  
[117] codetools_0.2-18              assertthat_0.2.1             
[119] rhdf5_2.36.0                  rjson_0.2.20                 
[121] withr_2.4.2                   GenomicAlignments_1.28.0     
[123] Rsamtools_2.8.0               GenomeInfoDbData_1.2.6       
[125] hms_1.1.0                     grid_4.1.0                   
[127] beachmat_2.8.0                rmarkdown_2.8                
[129] DelayedMatrixStats_1.14.0     Rtsne_0.15                   
[131] shiny_1.6.0                   ggbeeswarm_0.6.0             
[133] restfulr_0.0.13              
```
</div>