if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("SingleCellMultiModal")
library(SingleCellMultiModal)
library(MultiAssayExperiment)
library(scran)
library(scater)
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.
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 cache
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 files
where the cells are the same in both assays:
upsetSamples(mae)
Columns:
Lymphoid
or Myeloid
originThe 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.. Myeloid
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)