PBMCs profiled with the Chromium Single Cell Multiome ATAC + Gene Expression from 10x
6 May 2025
Source:vignettes/scMultiome.Rmd
scMultiome.RmdInstallation
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("SingleCellMultiModal")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.
Downloading datasets
The user can see the available dataset by using the default options
mae <- scMultiome("pbmc_10x", modes = "*", 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## Loading required namespace: HDF5Array## 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 cacheExploring 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)
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
## 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.. MyeloidRNA 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 10032
names(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)
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)
Suggested software for the downstream analysis
These are my personal recommendations of R-based analysis software:
sessionInfo
## R version 4.5.0 (2025-04-11)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 24.04.2 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.26.so; LAPACK version 3.12.0
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
## [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: Etc/UTC
## tzcode source: system (glibc)
##
## attached base packages:
## [1] stats4 stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] scater_1.36.0 ggplot2_3.5.2
## [3] scran_1.36.0 scuttle_1.18.0
## [5] SingleCellExperiment_1.30.0 SingleCellMultiModal_1.20.1
## [7] MultiAssayExperiment_1.34.0 SummarizedExperiment_1.38.1
## [9] Biobase_2.68.0 GenomicRanges_1.60.0
## [11] GenomeInfoDb_1.44.0 IRanges_2.42.0
## [13] S4Vectors_0.46.0 BiocGenerics_0.54.0
## [15] generics_0.1.3 MatrixGenerics_1.20.0
## [17] matrixStats_1.5.0 BiocStyle_2.36.0
##
## loaded via a namespace (and not attached):
## [1] RColorBrewer_1.1-3 jsonlite_2.0.0 magrittr_2.0.3
## [4] ggbeeswarm_0.7.2 magick_2.8.6 farver_2.1.2
## [7] rmarkdown_2.29 fs_1.6.6 ragg_1.4.0
## [10] vctrs_0.6.5 memoise_2.0.1 htmltools_0.5.8.1
## [13] S4Arrays_1.8.0 BiocBaseUtils_1.10.0 AnnotationHub_3.16.0
## [16] curl_6.2.2 BiocNeighbors_2.2.0 Rhdf5lib_1.30.0
## [19] SparseArray_1.8.0 rhdf5_2.52.0 sass_0.4.10
## [22] bslib_0.9.0 htmlwidgets_1.6.4 desc_1.4.3
## [25] plyr_1.8.9 cachem_1.1.0 igraph_2.1.4
## [28] mime_0.13 lifecycle_1.0.4 pkgconfig_2.0.3
## [31] rsvd_1.0.5 Matrix_1.7-3 R6_2.6.1
## [34] fastmap_1.2.0 GenomeInfoDbData_1.2.14 digest_0.6.37
## [37] AnnotationDbi_1.70.0 dqrng_0.4.1 irlba_2.3.5.1
## [40] ExperimentHub_2.16.0 textshaping_1.0.1 RSQLite_2.3.11
## [43] beachmat_2.24.0 filelock_1.0.3 labeling_0.4.3
## [46] httr_1.4.7 abind_1.4-8 compiler_4.5.0
## [49] bit64_4.6.0-1 withr_3.0.2 BiocParallel_1.42.0
## [52] viridis_0.6.5 DBI_1.2.3 UpSetR_1.4.0
## [55] HDF5Array_1.36.0 rappdirs_0.3.3 DelayedArray_0.34.1
## [58] rjson_0.2.23 bluster_1.18.0 tools_4.5.0
## [61] vipor_0.4.7 beeswarm_0.4.0 glue_1.8.0
## [64] h5mread_1.0.0 rhdf5filters_1.20.0 grid_4.5.0
## [67] cluster_2.1.8.1 gtable_0.3.6 BiocSingular_1.24.0
## [70] ScaledMatrix_1.16.0 metapod_1.16.0 XVector_0.48.0
## [73] RcppAnnoy_0.0.22 ggrepel_0.9.6 BiocVersion_3.21.1
## [76] pillar_1.10.2 limma_3.64.0 dplyr_1.1.4
## [79] BiocFileCache_2.16.0 lattice_0.22-7 bit_4.6.0
## [82] tidyselect_1.2.1 locfit_1.5-9.12 Biostrings_2.76.0
## [85] knitr_1.50 gridExtra_2.3 bookdown_0.43
## [88] edgeR_4.6.1 xfun_0.52 statmod_1.5.0
## [91] UCSC.utils_1.4.0 yaml_2.3.10 evaluate_1.0.3
## [94] codetools_0.2-20 tibble_3.2.1 BiocManager_1.30.25
## [97] cli_3.6.5 uwot_0.2.3 systemfonts_1.2.3
## [100] jquerylib_0.1.4 Rcpp_1.0.14 dbplyr_2.5.0
## [103] png_0.1-8 parallel_4.5.0 pkgdown_2.1.2
## [106] blob_1.2.4 SpatialExperiment_1.18.0 viridisLite_0.4.2
## [109] scales_1.4.0 purrr_1.0.4 crayon_1.5.3
## [112] rlang_1.1.6 KEGGREST_1.48.0