# 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 Section \@ref(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.0.4 (2021-02-15)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 20.04.2 LTS

Matrix products: default
BLAS:   /home/biocbuild/bbs-3.12-books/R/lib/libRblas.so
LAPACK: /home/biocbuild/bbs-3.12-books/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.18.5                EnsDb.Hsapiens.v86_2.99.0  
 [3] ensembldb_2.14.0            AnnotationFilter_1.14.0    
 [5] GenomicFeatures_1.42.2      AnnotationDbi_1.52.0       
 [7] scater_1.18.6               ggplot2_3.3.3              
 [9] DropletUtils_1.10.3         SingleCellExperiment_1.12.0
[11] SummarizedExperiment_1.20.0 Biobase_2.50.0             
[13] GenomicRanges_1.42.0        GenomeInfoDb_1.26.4        
[15] IRanges_2.24.1              S4Vectors_0.28.1           
[17] BiocGenerics_0.36.0         MatrixGenerics_1.2.1       
[19] matrixStats_0.58.0          DropletTestFiles_1.0.0     
[21] BiocStyle_2.18.1            rebook_1.0.0               

loaded via a namespace (and not attached):
  [1] AnnotationHub_2.22.0          BiocFileCache_1.14.0         
  [3] igraph_1.2.6                  lazyeval_0.2.2               
  [5] BiocParallel_1.24.1           digest_0.6.27                
  [7] htmltools_0.5.1.1             viridis_0.5.1                
  [9] fansi_0.4.2                   magrittr_2.0.1               
 [11] memoise_2.0.0                 limma_3.46.0                 
 [13] Biostrings_2.58.0             R.utils_2.10.1               
 [15] askpass_1.1                   prettyunits_1.1.1            
 [17] colorspace_2.0-0              blob_1.2.1                   
 [19] rappdirs_0.3.3                xfun_0.22                    
 [21] dplyr_1.0.5                   callr_3.5.1                  
 [23] crayon_1.4.1                  RCurl_1.98-1.3               
 [25] jsonlite_1.7.2                graph_1.68.0                 
 [27] glue_1.4.2                    gtable_0.3.0                 
 [29] zlibbioc_1.36.0               XVector_0.30.0               
 [31] DelayedArray_0.16.2           BiocSingular_1.6.0           
 [33] Rhdf5lib_1.12.1               HDF5Array_1.18.1             
 [35] scales_1.1.1                  DBI_1.1.1                    
 [37] edgeR_3.32.1                  Rcpp_1.0.6                   
 [39] viridisLite_0.3.0             xtable_1.8-4                 
 [41] progress_1.2.2                dqrng_0.2.1                  
 [43] bit_4.0.4                     rsvd_1.0.3                   
 [45] httr_1.4.2                    FNN_1.1.3                    
 [47] ellipsis_0.3.1                farver_2.1.0                 
 [49] pkgconfig_2.0.3               XML_3.99-0.6                 
 [51] R.methodsS3_1.8.1             scuttle_1.0.4                
 [53] uwot_0.1.10                   CodeDepends_0.6.5            
 [55] sass_0.3.1                    dbplyr_2.1.0                 
 [57] locfit_1.5-9.4                utf8_1.2.1                   
 [59] labeling_0.4.2                tidyselect_1.1.0             
 [61] rlang_0.4.10                  later_1.1.0.1                
 [63] munsell_0.5.0                 BiocVersion_3.12.0           
 [65] tools_4.0.4                   cachem_1.0.4                 
 [67] generics_0.1.0                RSQLite_2.2.4                
 [69] ExperimentHub_1.16.0          evaluate_0.14                
 [71] stringr_1.4.0                 fastmap_1.1.0                
 [73] yaml_2.2.1                    processx_3.4.5               
 [75] knitr_1.31                    bit64_4.0.5                  
 [77] purrr_0.3.4                   sparseMatrixStats_1.2.1      
 [79] mime_0.10                     R.oo_1.24.0                  
 [81] xml2_1.3.2                    biomaRt_2.46.3               
 [83] compiler_4.0.4                beeswarm_0.3.1               
 [85] curl_4.3                      interactiveDisplayBase_1.28.0
 [87] statmod_1.4.35                tibble_3.1.0                 
 [89] bslib_0.2.4                   stringi_1.5.3                
 [91] highr_0.8                     ps_1.6.0                     
 [93] RSpectra_0.16-0               bluster_1.0.0                
 [95] lattice_0.20-41               ProtGenerics_1.22.0          
 [97] Matrix_1.3-2                  vctrs_0.3.6                  
 [99] pillar_1.5.1                  lifecycle_1.0.0              
[101] rhdf5filters_1.2.0            BiocManager_1.30.10          
[103] jquerylib_0.1.3               BiocNeighbors_1.8.2          
[105] cowplot_1.1.1                 bitops_1.0-6                 
[107] irlba_2.3.3                   httpuv_1.5.5                 
[109] rtracklayer_1.50.0            R6_2.5.0                     
[111] bookdown_0.21                 promises_1.2.0.1             
[113] gridExtra_2.3                 vipor_0.4.5                  
[115] codetools_0.2-18              assertthat_0.2.1             
[117] rhdf5_2.34.0                  openssl_1.4.3                
[119] withr_2.4.1                   GenomicAlignments_1.26.0     
[121] Rsamtools_2.6.0               GenomeInfoDbData_1.2.4       
[123] hms_1.0.0                     grid_4.0.4                   
[125] beachmat_2.6.4                rmarkdown_2.7                
[127] DelayedMatrixStats_1.12.3     Rtsne_0.15                   
[129] shiny_1.6.0                   ggbeeswarm_0.6.0             
```
</div>