PBMCs profiled with the Chromium Single Cell Multiome ATAC + Gene Expression from 10x
5 September 2023
1 Installation
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("SingleCellMultiModal")1.1 Load
library(SingleCellMultiModal)
library(MultiAssayExperiment)
library(scran)
library(scater)2 Description
This data set consists of about 10K Peripheral Blood Mononuclear Cells (PBMCs) derived from a single healthy donor. It is available from the 10x Genomics website.
Provided are the RNA expression counts quantified at the gene level and the
chromatin accessibility levels quantified at the peak level. Here we provide
the default peaks called by the CellRanger software. If you want to explore
other peak definitions or chromatin accessibility quantifications (at the
promoter level, etc.), you have download the fragments.tsv.gz file from the
10x Genomics website.
3 Downloading datasets
The user can see the available dataset by using the default options
mae <- scMultiome("pbmc_10x", mode = "*", dry.run = FALSE, format = "MTX")## Working on: pbmc_atac_se.rds## Working on: pbmc_atac.mtx.gz## Working on: pbmc_rna_se.rds## Working on: pbmc_rna.mtx.gz## Working on: pbmc_atac,
## pbmc_rna## see ?SingleCellMultiModal and browseVignettes('SingleCellMultiModal') for documentation## loading from cache## see ?SingleCellMultiModal and browseVignettes('SingleCellMultiModal') for documentation## loading from cache## Working on: pbmc_atac,
## pbmc_rna## see ?SingleCellMultiModal and browseVignettes('SingleCellMultiModal') for documentation## loading from cache## see ?SingleCellMultiModal and browseVignettes('SingleCellMultiModal') for documentation## loading from cache## Working on: pbmc_colData## Working on: pbmc_sampleMap## see ?SingleCellMultiModal and browseVignettes('SingleCellMultiModal') for documentation## loading from cache## see ?SingleCellMultiModal and browseVignettes('SingleCellMultiModal') for documentation## loading from cache4 Exploring the data structure
There are two assays: rna and atac, stored as
SingleCellExperiment
objects
mae## A MultiAssayExperiment object of 2 listed
## experiments with user-defined names and respective classes.
## Containing an ExperimentList class object of length 2:
## [1] atac: SingleCellExperiment with 108344 rows and 10032 columns
## [2] rna: SingleCellExperiment with 36549 rows and 10032 columns
## Functionality:
## experiments() - obtain the ExperimentList instance
## colData() - the primary/phenotype DataFrame
## sampleMap() - the sample coordination DataFrame
## `$`, `[`, `[[` - extract colData columns, subset, or experiment
## *Format() - convert into a long or wide DataFrame
## assays() - convert ExperimentList to a SimpleList of matrices
## exportClass() - save data to flat fileswhere the cells are the same in both assays:
upsetSamples(mae)
4.1 Cell metadata
Columns:
- nCount_RNA: number of read counts
- nFeature_RNA: number of genes with at least one read count
- nCount_ATAC: number of ATAC read counts
- nFeature_ATAC: number of ATAC peaks with at least one read count
- celltype: The cell types have been annotated by the 10x Genomics R&D team using gene markers. They provide a rough characterisation of the cell type diversity, but keep in mind that they are not ground truth labels.
- broad_celltype:
LymphoidorMyeloidorigin
The cells have not been QC-ed, choosing a minimum number of genes/peaks per cell depends is left to you! In addition, there are further quality control criteria that you may want to apply, including mitochondrial coverage, fraction of reads overlapping ENCODE Blacklisted regions, Transcription start site enrichment, etc. See suggestions below for software that can perform a semi-automated quality control pipeline
head(colData(mae))## DataFrame with 6 rows and 6 columns
## nCount_RNA nFeature_RNA nCount_ATAC nFeature_ATAC
## <integer> <integer> <integer> <integer>
## AAACAGCCAAGGAATC 8380 3308 55582 13878
## AAACAGCCAATCCCTT 3771 1896 20495 7253
## AAACAGCCAATGCGCT 6876 2904 16674 6528
## AAACAGCCAGTAGGTG 7614 3061 39454 11633
## AAACAGCCAGTTTACG 3633 1691 20523 7245
## AAACAGCCATCCAGGT 7782 3028 22412 8602
## celltype broad_celltype
## <character> <character>
## AAACAGCCAAGGAATC naive CD4 T cells Lymphoid
## AAACAGCCAATCCCTT memory CD4 T cells Lymphoid
## AAACAGCCAATGCGCT naive CD4 T cells Lymphoid
## AAACAGCCAGTAGGTG naive CD4 T cells Lymphoid
## AAACAGCCAGTTTACG memory CD4 T cells Lymphoid
## AAACAGCCATCCAGGT non-classical monocy.. Myeloid4.2 RNA expression
The RNA expression consists of 36,549 genes and 10,032 cells, stored using
the dgCMatrix sparse matrix format
dim(experiments(mae)[["rna"]])## [1] 36549 10032names(experiments(mae))## [1] "atac" "rna"Let’s do some standard dimensionality reduction plot:
sce.rna <- experiments(mae)[["rna"]]
# Normalisation
sce.rna <- logNormCounts(sce.rna)
# Feature selection
decomp <- modelGeneVar(sce.rna)
hvgs <- rownames(decomp)[decomp$mean>0.01 & decomp$p.value <= 0.05]
sce.rna <- sce.rna[hvgs,]
# PCA
sce.rna <- runPCA(sce.rna, ncomponents = 25)
# UMAP
set.seed(42)
sce.rna <- runUMAP(sce.rna, dimred="PCA", n_neighbors = 25, min_dist = 0.3)
plotUMAP(sce.rna, colour_by="celltype", point_size=0.5, point_alpha=1)
4.3 Chromatin Accessibility
The ATAC expression consists of 108,344 peaks and 10,032 cells:
dim(experiments(mae)[["atac"]])## [1] 108344 10032Let’s do some standard dimensionality reduction plot. Note that scATAC-seq data is sparser than scRNA-seq, almost binary. The log normalisation + PCA approach that scater implements for scRNA-seq is not a good strategy for scATAC-seq data. Topic modelling or TFIDF+SVD are a better strategy. Please see the package recommendations below.
sce.atac <- experiments(mae)[["atac"]]
# Normalisation
sce.atac <- logNormCounts(sce.atac)
# Feature selection
decomp <- modelGeneVar(sce.atac)
hvgs <- rownames(decomp)[decomp$mean>0.25]
sce.atac <- sce.atac[hvgs,]
# PCA
sce.atac <- runPCA(sce.atac, ncomponents = 25)
# UMAP
set.seed(42)
sce.atac <- runUMAP(sce.atac, dimred="PCA", n_neighbors = 25, min_dist = 0.3)
plotUMAP(sce.atac, colour_by="celltype", point_size=0.5, point_alpha=1)
5 Suggested software for the downstream analysis
These are my personal recommendations of R-based analysis software:
6 sessionInfo
sessionInfo()## 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.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] scater_1.28.0 ggplot2_3.4.3
## [3] scran_1.28.2 scuttle_1.10.2
## [5] rhdf5_2.44.0 SingleCellExperiment_1.22.0
## [7] RaggedExperiment_1.24.2 SingleCellMultiModal_1.12.3
## [9] MultiAssayExperiment_1.26.0 SummarizedExperiment_1.30.2
## [11] Biobase_2.60.0 GenomicRanges_1.52.0
## [13] GenomeInfoDb_1.36.2 IRanges_2.34.1
## [15] S4Vectors_0.38.1 BiocGenerics_0.46.0
## [17] MatrixGenerics_1.12.3 matrixStats_1.0.0
## [19] BiocStyle_2.28.0
##
## loaded via a namespace (and not attached):
## [1] jsonlite_1.8.7 magrittr_2.0.3
## [3] ggbeeswarm_0.7.2 magick_2.7.5
## [5] farver_2.1.1 rmarkdown_2.24
## [7] zlibbioc_1.46.0 vctrs_0.6.3
## [9] memoise_2.0.1 DelayedMatrixStats_1.22.6
## [11] RCurl_1.98-1.12 htmltools_0.5.6
## [13] S4Arrays_1.0.6 BiocBaseUtils_1.2.0
## [15] AnnotationHub_3.8.0 curl_5.0.2
## [17] BiocNeighbors_1.18.0 Rhdf5lib_1.22.0
## [19] sass_0.4.7 bslib_0.5.1
## [21] plyr_1.8.8 cachem_1.0.8
## [23] igraph_1.5.1 mime_0.12
## [25] lifecycle_1.0.3 pkgconfig_2.0.3
## [27] rsvd_1.0.5 Matrix_1.6-1
## [29] R6_2.5.1 fastmap_1.1.1
## [31] GenomeInfoDbData_1.2.10 shiny_1.7.5
## [33] digest_0.6.33 colorspace_2.1-0
## [35] AnnotationDbi_1.62.2 dqrng_0.3.1
## [37] irlba_2.3.5.1 ExperimentHub_2.8.1
## [39] RSQLite_2.3.1 beachmat_2.16.0
## [41] labeling_0.4.3 filelock_1.0.2
## [43] fansi_1.0.4 httr_1.4.7
## [45] abind_1.4-5 compiler_4.3.1
## [47] bit64_4.0.5 withr_2.5.0
## [49] BiocParallel_1.34.2 viridis_0.6.4
## [51] DBI_1.1.3 UpSetR_1.4.0
## [53] highr_0.10 HDF5Array_1.28.1
## [55] R.utils_2.12.2 rappdirs_0.3.3
## [57] DelayedArray_0.26.7 rjson_0.2.21
## [59] bluster_1.10.0 tools_4.3.1
## [61] vipor_0.4.5 beeswarm_0.4.0
## [63] interactiveDisplayBase_1.38.0 httpuv_1.6.11
## [65] R.oo_1.25.0 glue_1.6.2
## [67] rhdf5filters_1.12.1 promises_1.2.1
## [69] grid_4.3.1 cluster_2.1.4
## [71] generics_0.1.3 gtable_0.3.4
## [73] R.methodsS3_1.8.2 BiocSingular_1.16.0
## [75] ScaledMatrix_1.8.1 metapod_1.8.0
## [77] utf8_1.2.3 XVector_0.40.0
## [79] RcppAnnoy_0.0.21 ggrepel_0.9.3
## [81] BiocVersion_3.17.1 pillar_1.9.0
## [83] limma_3.56.2 later_1.3.1
## [85] dplyr_1.1.3 BiocFileCache_2.8.0
## [87] lattice_0.21-8 bit_4.0.5
## [89] tidyselect_1.2.0 locfit_1.5-9.8
## [91] Biostrings_2.68.1 knitr_1.43
## [93] gridExtra_2.3 bookdown_0.35
## [95] edgeR_3.42.4 xfun_0.40
## [97] statmod_1.5.0 DropletUtils_1.20.0
## [99] yaml_2.3.7 evaluate_0.21
## [101] codetools_0.2-19 tibble_3.2.1
## [103] BiocManager_1.30.22 cli_3.6.1
## [105] uwot_0.1.16 xtable_1.8-4
## [107] munsell_0.5.0 jquerylib_0.1.4
## [109] Rcpp_1.0.11 dbplyr_2.3.3
## [111] png_0.1-8 parallel_4.3.1
## [113] ellipsis_0.3.2 blob_1.2.4
## [115] sparseMatrixStats_1.12.2 bitops_1.0-7
## [117] SpatialExperiment_1.10.0 viridisLite_0.4.2
## [119] scales_1.2.1 purrr_1.0.2
## [121] crayon_1.5.2 rlang_1.1.1
## [123] cowplot_1.1.1 KEGGREST_1.40.0
## [125] formatR_1.14