The estimated time for this lab is around 1h10.
ATAC-seq data may be obtained in isolation using a single-cell ATAC-seq protocol (e.g. 10X scATACseq) or in combination with gene expression data using a single-cell multiome protocole (e.g. 10X multiome and SHARE-seq).
Several packages are currently avaialble to process scATAC-seq data in R. These include Signac and ArchR. This lab will closely follow the processing steps outlined in Signac, which interfaces well with Seurat for single-cell analysis.
In this lab, we will process a PBMC single-cell ATAC-seq (scATAC-seq) dataset and perform preliminary analysis to assess quality of these data. The data for this lab comes from 10X Genomics. The dataset contains roughly ~ 10,000 cells.
Pre-downloaded data files are available in the ~/Share/scATAC/ folder. Notice how the count matrix is in a .h5 format. We have already encountered this format in Lab3. Back then, we imported it with DropletUtils::read10xCounts.
Does this function work here?
R
scatac <- DropletUtils::read10xCounts("~/Share/scATAC/atac_v1_pbmc_10k_filtered_peak_bc_matrix.h5")
scatacclass: SingleCellExperiment
dim: 89796 8728
metadata(1): Samples
assays(1): counts
rownames(89796): chr1:565107-565550 chr1:569174-569639 ... chrY:59001782-59002175 chrY:59017143-59017246
rowData names(3): ID Symbol Type
colnames: NULL
colData names(2): Sample Barcode
reducedDimNames(0):
mainExpName: NULL
altExpNames(0):This works because 10X Genomics make sure to distribute files in .h5 format that are consistent accross single-cell sequencing methodologies.
However, the SingleCellExperiment obtained with this approach is not the most convenient, as it cannot natively leverage fragments file (see below).
Instead, we can create a Signac object, a flavour of Seurat objects.
Import counts matrix and feature annotations using an import function provided by Seurat.
Loading required package: SeuratObjectLoading required package: sp
Attaching package: 'SeuratObject'The following objects are masked from 'package:base':
intersect, tLoading required package: GenomicRangesLoading required package: stats4Loading required package: BiocGenericsLoading required package: generics
Attaching package: 'generics'The following objects are masked from 'package:base':
as.difftime, as.factor, as.ordered, intersect, is.element, setdiff, setequal, union
Attaching package: 'BiocGenerics'The following objects are masked from 'package:stats':
IQR, mad, sd, var, xtabsThe following objects are masked from 'package:base':
anyDuplicated, aperm, append, as.data.frame, basename, cbind, colnames, dirname, do.call, duplicated, eval, evalq, Filter, Find, get, grep, grepl, is.unsorted, lapply, Map, mapply, match, mget, order, paste,
pmax, pmax.int, pmin, pmin.int, Position, rank, rbind, Reduce, rownames, sapply, saveRDS, table, tapply, unique, unsplit, which.max, which.minLoading required package: S4Vectors
Attaching package: 'S4Vectors'The following object is masked from 'package:utils':
findMatchesThe following objects are masked from 'package:base':
expand.grid, I, unnameLoading required package: IRanges
Attaching package: 'IRanges'The following object is masked from 'package:sp':
%over%Loading required package: Seqinfolibrary(stringr)
cnts <- Read10X_h5('~/Share/scATAC/atac_v1_pbmc_10k_filtered_peak_bc_matrix.h5')
features <- import('~/Share/scATAC/atac_v1_pbmc_10k_peaks.bed')
features$peak <- as.character(features) |> str_replace(':', '-')
metadata <- read.csv(
file = "~/Share/scATAC/atac_v1_pbmc_10k_singlecell.csv",
header = TRUE,
row.names = 1
)How many accessible genomic segments were found in this dataset?
R
featuresGRanges object with 89796 ranges and 1 metadata column:
seqnames ranges strand | peak
<Rle> <IRanges> <Rle> | <character>
[1] chr1 565108-565550 * | chr1-565108-565550
[2] chr1 569175-569639 * | chr1-569175-569639
[3] chr1 713461-714823 * | chr1-713461-714823
[4] chr1 752423-753038 * | chr1-752423-753038
[5] chr1 762107-763359 * | chr1-762107-763359
... ... ... ... . ...
[89792] chrY 58993393-58993760 * | chrY-58993393-58993760
[89793] chrY 58994572-58994823 * | chrY-58994572-58994823
[89794] chrY 58996353-58997331 * | chrY-58996353-58997331
[89795] chrY 59001783-59002175 * | chrY-59001783-59002175
[89796] chrY 59017144-59017246 * | chrY-59017144-59017246
-------
seqinfo: 24 sequences from an unspecified genome; no seqlengthslength(features)[1] 89796Seurat objectThe next step is to aggregate counts and features into a ChromatinAssay, a scATAC-seq flavour of Seurat standard Assays. The documentation for ?CreateChromatinAssay indicates that the user can provide:
fragments file, corresponding to the full list of all unique fragments mapped across all single cells.ChromatinAssay, corresponding to gene annotations, promoter positions, etc. Such annotations can be generated from Ensembl.Generate human annotations from Ensembl using a parsing function from Seurat.
R
## - Get human gene annotations (hg19/GRCh37) to feed it into the future `ChromatinAssay`
BiocManager::install('EnsDb.Hsapiens.v75')'getOption("repos")' replaces Bioconductor standard repositories, see 'help("repositories", package = "BiocManager")' for details.
Replacement repositories:
CRAN: https://cloud.r-project.orgBioconductor version 3.22 (BiocManager 1.30.27), R 4.5.1 (2025-06-13)Warning: package(s) not installed when version(s) same as or greater than current; use `force = TRUE` to re-install: 'EnsDb.Hsapiens.v75'Installation paths not writeable, unable to update packages
path: /opt/R/4.5.1/lib/R/library
packages:
boot, mgcvOld packages: 'affy', 'affyio', 'alphashape3d', 'AnnotationForge', 'AnnotationHubData', 'apeglm', 'aplot', 'arrow', 'aws.s3', 'AzureStor', 'basilisk.utils', 'BiocBook', 'BiocCheck', 'BiocPkgTools', 'biocthis', 'biocViews', 'biomformat',
'broom.helpers', 'BrowserViz', 'bs4Dash', 'bumphunter', 'BumpyMatrix', 'cards', 'CGHbase', 'CGHcall', 'ChIPpeakAnno', 'ChIPseeker', 'clusterProfiler', 'collections', 'ComplexHeatmap', 'corrr', 'covr', 'crisprBase', 'crisprBowtie',
'crisprScore', 'crisprScoreData', 'csaw', 'curatedMetagenomicData', 'curatedTCGAData', 'data.tree', 'dbscan', 'dearseq', 'DECIPHER', 'decontam', 'decoupleR', 'densvis', 'derfinder', 'derfinderHelper', 'DESeq2', 'DFplyr', 'dials',
'diffHic', 'DNAcopy', 'DNAZooData', 'DOSE', 'DropletTestFiles', 'easylift', 'EBImage', 'ellmer', 'emmeans', 'EnrichedHeatmap', 'enrichplot', 'eulerr', 'excluderanges', 'ExperimentHubData', 'fgsea', 'fission', 'fourDNData', 'fresh',
'genefilter', 'GenomeInfoDbData', 'GenomicFiles', 'GenomicInteractions', 'GEOquery', 'GGally', 'ggbio', 'ggExtra', 'ggimage', 'ggplotify', 'ggpubr', 'ggsci', 'ggstar', 'ggstats', 'ggtangle', 'ggtree', 'ggtreeExtra', 'gld', 'glmmTMB',
'GO.db', 'GOSemSim', 'GOTHiC', 'gprofiler2', 'gsl', 'GSVA', 'gtrellis', 'GUniFrac', 'Gviz', 'hardhat', 'HCAData', 'HiCcompare', 'HiCDataLymphoblast', 'HiCDOC', 'HiCExperiment', 'HiContacts', 'HiContactsData', 'igvR', 'impute',
'InteractionSet', 'interactiveDisplayBase', 'labelled', 'latticeExtra', 'lava', 'leaflet', 'lemon', 'LimROTS', 'litedown', 'lobstr', 'macrophage', 'magick', 'mariner', 'marray', 'MAST', 'Matrix', 'mclust', 'megadepth', 'memes',
'MetaboCoreUtils', 'metagene2', 'mia', 'miaViz', 'microbiomeDataSets', 'MicrobiotaProcess', 'mixOmics', 'MOFA2', 'monaLisa', 'MotifDb', 'motifmatchr', 'msa', 'MsBackendMetaboLights', 'msdata', 'msmsEDA', 'msmsTests', 'MSnbase', 'MSnID',
'MultiAssayExperiment', 'multiHiCcompare', 'muscat', 'mzID', 'mzR', 'NOISeq', 'nucleR', 'nullranges', 'oct4', 'ontoProc', 'ORFik', 'org.Hs.eg.db', 'org.Mm.eg.db', 'org.Sc.sgd.db', 'OrganismDbi', 'Orthology.eg.db', 'parsnip',
'pcaMethods', 'pdftools', 'periodicDNA', 'philentropy', 'philr', 'phyloseq', 'PKI', 'plotgardener', 'plyinteractions', 'plyranges', 'pracma', 'preprocessCore', 'PSMatch', 'QDNAseq', 'QFeatures', 'quarto', 'qvalue', 'RaggedExperiment',
'RBGL', 'Rbowtie', 'RcppParallel', 'rcrossref', 'REDseq', 'reformulas', 'regioneR', 'remaCor', 'Rfast', 'Rgraphviz', 'Rmpfr', 'ROC', 'rpx', 'rsconnect', 'rspm', 'rstatix', 'Rsubread', 'scatterpie', 'scmap', 'sechm', 'seriation',
'shadowtext', 'shinydashboardPlus', 'ShortRead', 'singscore', 'snowflakeauth', 'SpatialExperiment', 'Spectra', 'spelling', 'SRAdb', 'styler', 'subSeq', 'survey', 'svglite', 'svUnit', 'TCC', 'terra', 'tidybulk', 'tidyCoverage',
'tidySingleCellExperiment', 'tidySummarizedExperiment', 'timeDate', 'TMB', 'treeio', 'TreeSummarizedExperiment', 'tune', 'TxDb.Hsapiens.UCSC.hg19.knownGene', 'txdbmaker', 'tximeta', 'tximport', 'units', 'universalmotif',
'variancePartition', 'vegan', 'VplotR', 'vsn', 'waiter', 'workflows', 'Wrench', 'yulab.utils', 'zinbwave'annotations <- GetGRangesFromEnsDb(ensdb = EnsDb.Hsapiens.v75::EnsDb.Hsapiens.v75)Warning in .merge_two_Seqinfo_objects(x, y): The 2 combined objects have no sequence levels in common. (Use
suppressWarnings() to suppress this warning.)Warning in .merge_two_Seqinfo_objects(x, y): The 2 combined objects have no sequence levels in common. (Use
suppressWarnings() to suppress this warning.)
Warning in .merge_two_Seqinfo_objects(x, y): The 2 combined objects have no sequence levels in common. (Use
suppressWarnings() to suppress this warning.)
Warning in .merge_two_Seqinfo_objects(x, y): The 2 combined objects have no sequence levels in common. (Use
suppressWarnings() to suppress this warning.)
Warning in .merge_two_Seqinfo_objects(x, y): The 2 combined objects have no sequence levels in common. (Use
suppressWarnings() to suppress this warning.)
Warning in .merge_two_Seqinfo_objects(x, y): The 2 combined objects have no sequence levels in common. (Use
suppressWarnings() to suppress this warning.)
Warning in .merge_two_Seqinfo_objects(x, y): The 2 combined objects have no sequence levels in common. (Use
suppressWarnings() to suppress this warning.)
Warning in .merge_two_Seqinfo_objects(x, y): The 2 combined objects have no sequence levels in common. (Use
suppressWarnings() to suppress this warning.)
Warning in .merge_two_Seqinfo_objects(x, y): The 2 combined objects have no sequence levels in common. (Use
suppressWarnings() to suppress this warning.)
Warning in .merge_two_Seqinfo_objects(x, y): The 2 combined objects have no sequence levels in common. (Use
suppressWarnings() to suppress this warning.)
Warning in .merge_two_Seqinfo_objects(x, y): The 2 combined objects have no sequence levels in common. (Use
suppressWarnings() to suppress this warning.)
Warning in .merge_two_Seqinfo_objects(x, y): The 2 combined objects have no sequence levels in common. (Use
suppressWarnings() to suppress this warning.)
Warning in .merge_two_Seqinfo_objects(x, y): The 2 combined objects have no sequence levels in common. (Use
suppressWarnings() to suppress this warning.)
Warning in .merge_two_Seqinfo_objects(x, y): The 2 combined objects have no sequence levels in common. (Use
suppressWarnings() to suppress this warning.)
Warning in .merge_two_Seqinfo_objects(x, y): The 2 combined objects have no sequence levels in common. (Use
suppressWarnings() to suppress this warning.)
Warning in .merge_two_Seqinfo_objects(x, y): The 2 combined objects have no sequence levels in common. (Use
suppressWarnings() to suppress this warning.)
Warning in .merge_two_Seqinfo_objects(x, y): The 2 combined objects have no sequence levels in common. (Use
suppressWarnings() to suppress this warning.)
Warning in .merge_two_Seqinfo_objects(x, y): The 2 combined objects have no sequence levels in common. (Use
suppressWarnings() to suppress this warning.)
Warning in .merge_two_Seqinfo_objects(x, y): The 2 combined objects have no sequence levels in common. (Use
suppressWarnings() to suppress this warning.)
Warning in .merge_two_Seqinfo_objects(x, y): The 2 combined objects have no sequence levels in common. (Use
suppressWarnings() to suppress this warning.)
Warning in .merge_two_Seqinfo_objects(x, y): The 2 combined objects have no sequence levels in common. (Use
suppressWarnings() to suppress this warning.)
Warning in .merge_two_Seqinfo_objects(x, y): The 2 combined objects have no sequence levels in common. (Use
suppressWarnings() to suppress this warning.)
Warning in .merge_two_Seqinfo_objects(x, y): The 2 combined objects have no sequence levels in common. (Use
suppressWarnings() to suppress this warning.)
Warning in .merge_two_Seqinfo_objects(x, y): The 2 combined objects have no sequence levels in common. (Use
suppressWarnings() to suppress this warning.)GenomeInfoDb::seqlevelsStyle(annotations) <- 'UCSC'Create a ChromatinAssay using counts, features, fragments and annotations.
R
## - Create Chromatin Assay
assay <- Signac::CreateChromatinAssay(
counts = cnts,
ranges = features,
fragments = "~/Share/scATAC/atac_v1_pbmc_10k_fragments.tsv.gz",
annotation = annotations,
genome = "hg19",
min.cells = 10,
min.features = 10
)Computing hashassayChromatinAssay data with 87561 features for 8728 cells
Variable features: 0
Genome: hg19
Annotation present: TRUE
Motifs present: FALSE
Fragment files: 1 What are the dimensions of this object? Are they comparable to the count matrix? Comment.
It’s finally time to wrap the ChromatinAssay into a Seurat standard object! This is done using the CreateSeuratObject, as already covered in Lab6
Create a Seurat object from the ChromatinAssay object.
R
## - Create Seurat object
PBMC <- Seurat::CreateSeuratObject(
counts = assay,
assay = 'ATAC',
meta.data = metadata[metadata$is__cell_barcode == 1, ]
)
PBMCAn object of class Seurat
87561 features across 8728 samples within 1 assay
Active assay: ATAC (87561 features, 0 variable features)
2 layers present: counts, dataPBMC[['ATAC']]ChromatinAssay data with 87561 features for 8728 cells
Variable features: 0
Genome: hg19
Annotation present: TRUE
Motifs present: FALSE
Fragment files: 1 granges(PBMC)GRanges object with 87561 ranges and 1 metadata column:
seqnames ranges strand | peak
<Rle> <IRanges> <Rle> | <character>
[1] chr1 565108-565550 * | chr1-565108-565550
[2] chr1 569175-569639 * | chr1-569175-569639
[3] chr1 713461-714823 * | chr1-713461-714823
[4] chr1 752423-753038 * | chr1-752423-753038
[5] chr1 762107-763359 * | chr1-762107-763359
... ... ... ... . ...
[87557] chrY 58993393-58993760 * | chrY-58993393-58993760
[87558] chrY 58994572-58994823 * | chrY-58994572-58994823
[87559] chrY 58996353-58997331 * | chrY-58996353-58997331
[87560] chrY 59001783-59002175 * | chrY-59001783-59002175
[87561] chrY 59017144-59017246 * | chrY-59017144-59017246
-------
seqinfo: 24 sequences from an unspecified genome; no seqlengthsAnnotation(PBMC)GRanges object with 3072120 ranges and 5 metadata columns:
seqnames ranges strand | tx_id gene_name gene_id gene_biotype type
<Rle> <IRanges> <Rle> | <character> <character> <character> <character> <factor>
ENSE00001489430 chrX 192989-193061 + | ENST00000399012 PLCXD1 ENSG00000182378 protein_coding exon
ENSE00001536003 chrX 192991-193061 + | ENST00000484611 PLCXD1 ENSG00000182378 protein_coding exon
ENSE00002160563 chrX 193020-193061 + | ENST00000430923 PLCXD1 ENSG00000182378 protein_coding exon
ENSE00001750899 chrX 197722-197788 + | ENST00000445062 PLCXD1 ENSG00000182378 protein_coding exon
ENSE00001489388 chrX 197859-198351 + | ENST00000381657 PLCXD1 ENSG00000182378 protein_coding exon
... ... ... ... . ... ... ... ... ...
ENST00000361739 chrMT 7586-8269 + | ENST00000361739 MT-CO2 ENSG00000198712 protein_coding cds
ENST00000361789 chrMT 14747-15887 + | ENST00000361789 MT-CYB ENSG00000198727 protein_coding cds
ENST00000361851 chrMT 8366-8572 + | ENST00000361851 MT-ATP8 ENSG00000228253 protein_coding cds
ENST00000361899 chrMT 8527-9207 + | ENST00000361899 MT-ATP6 ENSG00000198899 protein_coding cds
ENST00000362079 chrMT 9207-9990 + | ENST00000362079 MT-CO3 ENSG00000198938 protein_coding cds
-------
seqinfo: 25 sequences (1 circular) from hg19 genomeThe fraction of reads in peaks (FRiP) is a good indicator of how well each cell was handled during scATACseq processing.
R
PBMC$FRiP <- PBMC$peak_region_fragments / PBMC$passed_filters
PBMC$nCount_ATAC <- colSums(GetAssayData(PBMC, layer = "counts"))
PBMC$nFeature_ATAC <- colSums(GetAssayData(PBMC, layer = "counts") > 0)
quantile(PBMC$FRiP, seq(0, 1, 0.1)) 0% 10% 20% 30% 40% 50% 60% 70% 80% 90% 100%
0.02770289 0.55266699 0.60895480 0.64086129 0.66583010 0.68520756 0.70194678 0.71993709 0.73801044 0.76004660 0.90383142 0% 10% 20% 30% 40% 50% 60% 70% 80% 90% 100%
1710.0 5079.1 7619.8 9522.5 11300.8 13367.5 15754.6 19005.5 23502.6 31316.2 249799.0 0% 10% 20% 30% 40% 50% 60% 70% 80% 90% 100%
654.0 2308.4 3349.4 4108.1 4765.0 5447.0 6234.2 7100.9 8241.6 9983.4 30313.0 Which analysis are the fragments required for, exactly?
Could we still perform normalization/clustering/annotation without them? And motif enrichment analysis?
Since we do have the fragments file at hand, most of the QC steps are available (e.g. TSSEnrichment, NucleosomeSignal or fragment size distribution). Let’s go through them one by one.
R
# compute nucleosome signal score per cell
PBMC <- NucleosomeSignal(object = PBMC)
PBMCAn object of class Seurat
87561 features across 8728 samples within 1 assay
Active assay: ATAC (87561 features, 0 variable features)
2 layers present: counts, data# compute TSS enrichment score per cell
PBMC <- Signac::TSSEnrichment(object = PBMC, fast = FALSE)Extracting TSS positionsFinding + strand cut sitesFinding - strand cut sitesComputing mean insertion frequency in flanking regionsNormalizing TSS scorePBMCAn object of class Seurat
87561 features across 8728 samples within 1 assay
Active assay: ATAC (87561 features, 0 variable features)
2 layers present: counts, dataThe TSSPlot function from Signac can be used to plot the fragment count per peak ~ TSS enrichment.
R
Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
ℹ Please use `linewidth` instead.
ℹ The deprecated feature was likely used in the Signac package.
Please report the issue at <https://github.com/stuart-lab/signac/issues>.
PBMC$high.tss <- ifelse(PBMC$TSS.enrichment > 2.5, 'High', 'Low')
TSSPlot(PBMC, group.by = 'high.tss') + NoLegend()
The FragmentHistogram function from Signac can be used to plot the fragment size distribution in peaks with different nucleosome signals.
Group genomic loci by their nucleosome signal (> 4 or < 4) and plot the fragment size distribution.
R
PBMC$nucleosome_group <- ifelse(PBMC$nucleosome_signal > 4, 'NS > 4', 'NS < 4')
FragmentHistogram(object = PBMC, group.by = 'nucleosome_group')Warning: Removed 92 rows containing non-finite outside the scale range (`stat_bin()`).Warning: Removed 4 rows containing missing values or values outside the scale range (`geom_bar()`).
Just like standard scRNAseq data, scATACseq data needs to be filtered to remove low-quality cells and peaks.
Filter the Seurat object (cells and features) as deemed appropriate.
R
An object of class Seurat
87561 features across 8231 samples within 1 assay
Active assay: ATAC (87561 features, 0 variable features)
2 layers present: counts, data## - Remove peaks with low coverage
PBMC <- PBMC[rowSums(GetAssayData(PBMC, layer = "counts")) > 10, ]
PBMC <- PBMC[rowSums(GetAssayData(PBMC, layer = "counts") > 0) > 10, ]
PBMCAn object of class Seurat
87212 features across 8231 samples within 1 assay
Active assay: ATAC (87212 features, 0 variable features)
2 layers present: counts, dataNow that the dataset is filtered, we can proceed to data normalization, dimensionality reduction and clustering. Compared to scRNAseq, scATACseq data is much sparser and has a higher dimensionality. This has important implications for the choice of normalization and dimensionality reduction methods. scATACseq data are generally normalized using a TF-IDF normalization, which is implemented in the RunTFIDF function in Signac, and dimensionality reduction is performed using LSI (Latent Semantic Indexing) implemented in the RunSVD function.
Normalize then further reduce the dimensionality for visualization purposes.
R
## - Normalize data
PBMC <- RunTFIDF(PBMC) Performing TF-IDF normalization## - Reduce dimensionality
PBMC <- FindTopFeatures(PBMC, min.cutoff = 'q50')
PBMC <- RunSVD(PBMC) Running SVDScaling cell embeddings## - Embed in 2D UMAP
PBMC <- RunUMAP(object = PBMC, reduction = 'lsi', dims = 2:30)Warning: The default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metric
To use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation'
This message will be shown once per session20:24:04 UMAP embedding parameters a = 0.9922 b = 1.11220:24:04 Read 8231 rows and found 29 numeric columns20:24:04 Using Annoy for neighbor search, n_neighbors = 3020:24:04 Building Annoy index with metric = cosine, n_trees = 500% 10 20 30 40 50 60 70 80 90 100%[----|----|----|----|----|----|----|----|----|----|**************************************************|
20:24:05 Writing NN index file to temp file /data/.cache/R//RtmpaFe87f/file16a08c51315069
20:24:05 Searching Annoy index using 1 thread, search_k = 3000
20:24:07 Annoy recall = 100%
20:24:08 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
20:24:09 Initializing from normalized Laplacian + noise (using RSpectra)
20:24:09 Commencing optimization for 500 epochs, with 316686 positive edges
20:24:09 Using rng type: pcg
20:24:16 Optimization finished## - Cluster cells
PBMC <- FindNeighbors(object = PBMC, reduction = 'lsi', dims = 2:30)Computing nearest neighbor graph
Computing SNNPBMC <- FindClusters(object = PBMC, verbose = FALSE, algorithm = 3)
## - Visualize data
DimPlot(PBMC, label = TRUE) + NoLegend()
What can you observe in the UMAP projection of the dataset? Comment on the separation of some cell types into different spatially-resolved clusters.
Signac’s GeneActivity() function require scATACseq fragment information. Since we have them, we can estimate a gene activity score for each gene in the annotations.
R
gene.activities <- GeneActivity(PBMC)Extracting gene coordinatesWarning in SingleFeatureMatrix(fragment = fragments[[x]], features = features, : 13 features are on seqnames not present in the fragment file. These will be removed.Extracting reads overlapping genomic regionsWe can now save this new object as an Assay in the PBMC object and normalize it.
R
PBMC[['RNA']] <- CreateAssayObject(counts = gene.activities)Warning: Non-unique features (rownames) present in the input matrix, making unique# - Normalize the new RNA assay.
PBMC <- NormalizeData(object = PBMC, assay = 'RNA', normalization.method = 'LogNormalize', scale.factor = median(PBMC$nCount_RNA))
PBMCAn object of class Seurat
107222 features across 8231 samples within 2 assays
Active assay: ATAC (87212 features, 43610 variable features)
2 layers present: counts, data
1 other assay present: RNA
2 dimensional reductions calculated: lsi, umapOne can now visualize expression of individual markers across clusters to infer cluster identity.
Assess the expression of known markers across clusters, projecting the cells in a 2D UMAP.
R
PBMC$annotation <- dplyr::case_when(
Idents(PBMC) %in% c(0, 5, 9, 12) ~ "Monocytes",
Idents(PBMC) %in% c(1, 2, 3, 4, 7, 10) ~ "T-cells",
Idents(PBMC) %in% c(6) ~ "NK-cells",
Idents(PBMC) %in% c(8, 11) ~ "B-cells",
Idents(PBMC) %in% c(13, 14) ~ "DC"
)
DimPlot(PBMC, group.by = 'annotation', label = TRUE) + NoLegend()
Finally, we can perform differential peak accessibility analysis to identify peaks that are differentially accessible between clusters. This is done using the FindMarkers function in Signac, which is similar to FindAllMarkers in Seurat.
R
Idents(PBMC) <- PBMC$annotation
da_peaks <- FindMarkers(
object = PBMC,
ident.1 = "Monocytes",
ident.2 = "B-cells",
test.use = 'LR',
latent.vars = 'nCount_ATAC'
)
VlnPlot(
object = PBMC,
features = rownames(da_peaks)[1],
pt.size = 0.1,
idents = c("Monocytes","B-cells")
)─ Session info ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
setting value
version R version 4.5.1 (2025-06-13)
os Ubuntu 24.04.2 LTS
system x86_64, linux-gnu
ui X11
language (EN)
collate en_IN.utf8
ctype en_IN.utf8
tz Europe/Paris
date 2025-11-14
pandoc 3.1.3 @ /usr/bin/ (via rmarkdown)
quarto 1.5.55 @ /opt/quarto/bin/quarto
─ Packages ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
package * version date (UTC) lib source
abind 1.4-8 2024-09-12 [1] CRAN (R 4.5.1)
AnnotationDbi 1.72.0 2025-10-29 [1] Bioconductor 3.22 (R 4.5.1)
AnnotationFilter 1.34.0 2025-10-29 [1] Bioconductor 3.22 (R 4.5.1)
backports 1.5.0 2024-05-23 [1] CRAN (R 4.5.1)
base64enc 0.1-3 2015-07-28 [1] CRAN (R 4.5.1)
beachmat 2.26.0 2025-10-29 [1] Bioconductor 3.22 (R 4.5.1)
beeswarm 0.4.0 2021-06-01 [1] CRAN (R 4.5.1)
Biobase 2.70.0 2025-10-29 [1] Bioconductor 3.22 (R 4.5.1)
BiocGenerics * 0.56.0 2025-10-29 [1] Bioconductor 3.22 (R 4.5.1)
BiocIO 1.20.0 2025-10-29 [1] Bioconductor 3.22 (R 4.5.1)
BiocManager 1.30.27 2025-11-14 [1] CRAN (R 4.5.1)
BiocParallel 1.44.0 2025-10-29 [1] Bioconductor 3.22 (R 4.5.1)
Biostrings 2.78.0 2025-10-29 [1] Bioconductor 3.22 (R 4.5.1)
biovizBase 1.58.0 2025-10-29 [1] Bioconductor 3.22 (R 4.5.1)
bit 4.6.0 2025-03-06 [1] CRAN (R 4.5.0)
bit64 4.6.0-1 2025-01-16 [1] CRAN (R 4.5.0)
bitops 1.0-9 2024-10-03 [1] CRAN (R 4.5.1)
blob 1.2.4 2023-03-17 [1] CRAN (R 4.5.1)
BSgenome 1.78.0 2025-10-29 [1] Bioconductor 3.22 (R 4.5.1)
cachem 1.1.0 2024-05-16 [1] CRAN (R 4.5.1)
checkmate 2.3.3 2025-08-18 [1] CRAN (R 4.5.1)
cigarillo 1.0.0 2025-10-29 [1] Bioconductor 3.22 (R 4.5.1)
cli 3.6.5 2025-04-23 [1] CRAN (R 4.5.1)
cluster 2.1.8.1 2025-03-12 [2] CRAN (R 4.5.1)
codetools 0.2-20 2024-03-31 [2] CRAN (R 4.5.1)
colorspace 2.1-2 2025-09-22 [1] CRAN (R 4.5.1)
cowplot 1.2.0 2025-07-07 [1] CRAN (R 4.5.1)
crayon 1.5.3 2024-06-20 [1] CRAN (R 4.5.1)
curl 7.0.0 2025-08-19 [1] CRAN (R 4.5.1)
data.table 1.17.8 2025-07-10 [1] CRAN (R 4.5.1)
DBI 1.2.3 2024-06-02 [1] CRAN (R 4.5.1)
DelayedArray 0.36.0 2025-10-29 [1] Bioconductor 3.22 (R 4.5.1)
DelayedMatrixStats 1.32.0 2025-10-29 [1] Bioconductor 3.22 (R 4.5.1)
deldir 2.0-4 2024-02-28 [1] CRAN (R 4.5.1)
devtools 2.4.6 2025-10-03 [1] CRAN (R 4.5.1)
dichromat 2.0-0.1 2022-05-02 [1] CRAN (R 4.5.1)
digest 0.6.38 2025-11-12 [1] CRAN (R 4.5.1)
dotCall64 1.2 2024-10-04 [1] CRAN (R 4.5.1)
dplyr 1.1.4 2023-11-17 [1] CRAN (R 4.5.1)
dqrng 0.4.1 2024-05-28 [1] CRAN (R 4.5.1)
DropletUtils 1.30.0 2025-10-29 [1] Bioconductor 3.22 (R 4.5.1)
edgeR 4.8.0 2025-10-29 [1] Bioconductor 3.22 (R 4.5.1)
ellipsis 0.3.2 2021-04-29 [1] CRAN (R 4.5.1)
EnsDb.Hsapiens.v75 2.99.0 2025-07-16 [1] Bioconductor
ensembldb 2.34.0 2025-10-29 [1] Bioconductor 3.22 (R 4.5.1)
evaluate 1.0.5 2025-08-27 [1] CRAN (R 4.5.1)
farver 2.1.2 2024-05-13 [1] CRAN (R 4.5.1)
fastDummies 1.7.5 2025-01-20 [1] CRAN (R 4.5.0)
fastmap 1.2.0 2024-05-15 [1] CRAN (R 4.5.1)
fastmatch 1.1-6 2024-12-23 [1] CRAN (R 4.5.0)
fitdistrplus 1.2-4 2025-07-03 [1] CRAN (R 4.5.1)
foreign 0.8-90 2025-03-31 [2] CRAN (R 4.5.1)
Formula 1.2-5 2023-02-24 [1] CRAN (R 4.5.1)
fs 1.6.6 2025-04-12 [1] CRAN (R 4.5.1)
future * 1.67.0 2025-07-29 [1] CRAN (R 4.5.1)
future.apply 1.20.0 2025-06-06 [1] CRAN (R 4.5.1)
generics * 0.1.4 2025-05-09 [1] CRAN (R 4.5.1)
GenomeInfoDb 1.46.0 2025-10-29 [1] Bioconductor 3.22 (R 4.5.1)
GenomicAlignments 1.46.0 2025-10-29 [1] Bioconductor 3.22 (R 4.5.1)
GenomicFeatures 1.62.0 2025-10-29 [1] Bioconductor 3.22 (R 4.5.1)
GenomicRanges * 1.62.0 2025-10-29 [1] Bioconductor 3.22 (R 4.5.1)
ggbeeswarm 0.7.2 2023-04-29 [1] CRAN (R 4.5.1)
ggplot2 4.0.1 2025-11-14 [1] CRAN (R 4.5.1)
ggrastr 1.0.2 2023-06-01 [1] CRAN (R 4.5.1)
ggrepel 0.9.6 2024-09-07 [1] CRAN (R 4.5.1)
ggridges 0.5.7 2025-08-27 [1] CRAN (R 4.5.1)
globals 0.18.0 2025-05-08 [1] CRAN (R 4.5.1)
glue 1.8.0 2024-09-30 [1] CRAN (R 4.5.1)
goftest 1.2-3 2021-10-07 [1] CRAN (R 4.5.1)
gridExtra 2.3 2017-09-09 [1] CRAN (R 4.5.1)
gtable 0.3.6 2024-10-25 [1] CRAN (R 4.5.1)
h5mread 1.2.0 2025-10-29 [1] Bioconductor 3.22 (R 4.5.1)
HDF5Array 1.38.0 2025-10-29 [1] Bioconductor 3.22 (R 4.5.1)
hdf5r 1.3.12 2025-01-20 [1] CRAN (R 4.5.1)
Hmisc 5.2-4 2025-10-05 [1] CRAN (R 4.5.1)
htmlTable 2.4.3 2024-07-21 [1] CRAN (R 4.5.1)
htmltools 0.5.8.1 2024-04-04 [1] CRAN (R 4.5.1)
htmlwidgets 1.6.4 2023-12-06 [1] CRAN (R 4.5.1)
httpuv 1.6.16 2025-04-16 [1] CRAN (R 4.5.1)
httr 1.4.7 2023-08-15 [1] CRAN (R 4.5.1)
ica 1.0-3 2022-07-08 [1] CRAN (R 4.5.1)
igraph 2.2.1 2025-10-27 [1] CRAN (R 4.5.1)
IRanges * 2.44.0 2025-10-29 [1] Bioconductor 3.22 (R 4.5.1)
irlba 2.3.5.1 2022-10-03 [1] CRAN (R 4.5.1)
jsonlite 2.0.0 2025-03-27 [1] CRAN (R 4.5.1)
KEGGREST 1.50.0 2025-10-29 [1] Bioconductor 3.22 (R 4.5.1)
KernSmooth 2.23-26 2025-01-01 [2] CRAN (R 4.5.1)
knitr 1.50 2025-03-16 [1] CRAN (R 4.5.1)
labeling 0.4.3 2023-08-29 [1] CRAN (R 4.5.1)
later 1.4.4 2025-08-27 [1] CRAN (R 4.5.1)
lattice 0.22-7 2025-04-02 [2] CRAN (R 4.5.1)
lazyeval 0.2.2 2019-03-15 [1] CRAN (R 4.5.1)
lifecycle 1.0.4 2023-11-07 [1] CRAN (R 4.5.1)
limma 3.66.0 2025-10-29 [1] Bioconductor 3.22 (R 4.5.1)
listenv 0.10.0 2025-11-02 [1] CRAN (R 4.5.1)
lmtest 0.9-40 2022-03-21 [1] CRAN (R 4.5.1)
locfit 1.5-9.12 2025-03-05 [1] CRAN (R 4.5.0)
magrittr 2.0.4 2025-09-12 [1] CRAN (R 4.5.1)
MASS 7.3-65 2025-02-28 [2] CRAN (R 4.5.1)
Matrix 1.7-3 2025-03-11 [1] CRAN (R 4.5.1)
MatrixGenerics 1.22.0 2025-10-29 [1] Bioconductor 3.22 (R 4.5.1)
matrixStats 1.5.0 2025-01-07 [1] CRAN (R 4.5.0)
memoise 2.0.1 2021-11-26 [1] CRAN (R 4.5.1)
mime 0.13 2025-03-17 [1] CRAN (R 4.5.1)
miniUI 0.1.2 2025-04-17 [1] CRAN (R 4.5.1)
nlme 3.1-168 2025-03-31 [2] CRAN (R 4.5.1)
nnet 7.3-20 2025-01-01 [2] CRAN (R 4.5.1)
otel 0.2.0 2025-08-29 [1] CRAN (R 4.5.1)
parallelly 1.45.1 2025-07-24 [1] CRAN (R 4.5.1)
patchwork 1.3.2 2025-08-25 [1] CRAN (R 4.5.1)
pbapply 1.7-4 2025-07-20 [1] CRAN (R 4.5.1)
pillar 1.11.1 2025-09-17 [1] CRAN (R 4.5.1)
pkgbuild 1.4.8 2025-05-26 [1] CRAN (R 4.5.1)
pkgconfig 2.0.3 2019-09-22 [1] CRAN (R 4.5.1)
pkgload 1.4.1 2025-09-23 [1] CRAN (R 4.5.1)
plotly 4.11.0 2025-06-19 [1] CRAN (R 4.5.1)
plyr 1.8.9 2023-10-02 [1] CRAN (R 4.5.1)
png 0.1-8 2022-11-29 [1] CRAN (R 4.5.1)
polyclip 1.10-7 2024-07-23 [1] CRAN (R 4.5.1)
progressr 0.18.0 2025-11-06 [1] CRAN (R 4.5.1)
promises 1.5.0 2025-11-01 [1] CRAN (R 4.5.1)
ProtGenerics 1.42.0 2025-10-29 [1] Bioconductor 3.22 (R 4.5.1)
purrr 1.2.0 2025-11-04 [1] CRAN (R 4.5.1)
R.methodsS3 1.8.2 2022-06-13 [1] CRAN (R 4.5.1)
R.oo 1.27.1 2025-05-02 [1] CRAN (R 4.5.1)
R.utils 2.13.0 2025-02-24 [1] CRAN (R 4.5.0)
R6 2.6.1 2025-02-15 [1] CRAN (R 4.5.0)
RANN 2.6.2 2024-08-25 [1] CRAN (R 4.5.1)
RColorBrewer 1.1-3 2022-04-03 [1] CRAN (R 4.5.1)
Rcpp 1.1.0 2025-07-02 [1] CRAN (R 4.5.1)
RcppAnnoy 0.0.22 2024-01-23 [1] CRAN (R 4.5.1)
RcppHNSW 0.6.0 2024-02-04 [1] CRAN (R 4.5.1)
RcppRoll 0.3.1 2024-07-07 [1] CRAN (R 4.5.1)
RCurl 1.98-1.17 2025-03-22 [1] CRAN (R 4.5.1)
remotes 2.5.0 2024-03-17 [1] CRAN (R 4.5.1)
reshape2 1.4.5 2025-11-12 [1] CRAN (R 4.5.1)
restfulr 0.0.16 2025-06-27 [1] CRAN (R 4.5.1)
reticulate 1.44.0 2025-10-25 [1] CRAN (R 4.5.1)
rhdf5 2.54.0 2025-10-29 [1] Bioconductor 3.22 (R 4.5.1)
rhdf5filters 1.22.0 2025-10-29 [1] Bioconductor 3.22 (R 4.5.1)
Rhdf5lib 1.32.0 2025-10-29 [1] Bioconductor 3.22 (R 4.5.1)
rjson 0.2.23 2024-09-16 [1] CRAN (R 4.5.1)
rlang 1.1.6 2025-04-11 [1] CRAN (R 4.5.1)
rmarkdown 2.30 2025-09-28 [1] CRAN (R 4.5.1)
ROCR 1.0-11 2020-05-02 [1] CRAN (R 4.5.1)
rpart 4.1.24 2025-01-07 [2] CRAN (R 4.5.1)
Rsamtools 2.26.0 2025-10-29 [1] Bioconductor 3.22 (R 4.5.1)
RSpectra 0.16-2 2024-07-18 [1] CRAN (R 4.5.1)
RSQLite 2.4.4 2025-11-10 [1] CRAN (R 4.5.1)
rstudioapi 0.17.1 2024-10-22 [1] CRAN (R 4.5.1)
rtracklayer * 1.70.0 2025-10-29 [1] Bioconductor 3.22 (R 4.5.1)
Rtsne 0.17 2023-12-07 [1] CRAN (R 4.5.1)
S4Arrays 1.10.0 2025-10-29 [1] Bioconductor 3.22 (R 4.5.1)
S4Vectors * 0.48.0 2025-10-29 [1] Bioconductor 3.22 (R 4.5.1)
S7 0.2.0 2024-11-07 [1] CRAN (R 4.5.1)
scales 1.4.0 2025-04-24 [1] CRAN (R 4.5.1)
scattermore 1.2 2023-06-12 [1] CRAN (R 4.5.1)
sctransform 0.4.2 2025-04-30 [1] CRAN (R 4.5.1)
scuttle 1.20.0 2025-10-30 [1] Bioconductor 3.22 (R 4.5.1)
Seqinfo * 1.0.0 2025-10-29 [1] Bioconductor 3.22 (R 4.5.1)
sessioninfo 1.2.3 2025-02-05 [1] CRAN (R 4.5.0)
Seurat * 5.3.1 2025-10-29 [1] CRAN (R 4.5.1)
SeuratObject * 5.2.0 2025-08-27 [1] CRAN (R 4.5.1)
shiny 1.11.1 2025-07-03 [1] CRAN (R 4.5.1)
Signac * 1.16.0 2025-10-10 [1] CRAN (R 4.5.1)
SingleCellExperiment 1.32.0 2025-10-29 [1] Bioconductor 3.22 (R 4.5.1)
sp * 2.2-0 2025-02-01 [1] CRAN (R 4.5.0)
spam 2.11-1 2025-01-20 [1] CRAN (R 4.5.0)
SparseArray 1.10.1 2025-10-31 [1] Bioconductor 3.22 (R 4.5.1)
sparseMatrixStats 1.22.0 2025-10-29 [1] Bioconductor 3.22 (R 4.5.1)
spatstat.data 3.1-9 2025-10-18 [1] CRAN (R 4.5.1)
spatstat.explore 3.5-3 2025-09-22 [1] CRAN (R 4.5.1)
spatstat.geom 3.6-0 2025-09-20 [1] CRAN (R 4.5.1)
spatstat.random 3.4-2 2025-09-21 [1] CRAN (R 4.5.1)
spatstat.sparse 3.1-0 2024-06-21 [1] CRAN (R 4.5.1)
spatstat.univar 3.1-4 2025-07-13 [1] CRAN (R 4.5.1)
spatstat.utils 3.2-0 2025-09-20 [1] CRAN (R 4.5.1)
statmod 1.5.1 2025-10-09 [1] CRAN (R 4.5.1)
stringi 1.8.7 2025-03-27 [1] CRAN (R 4.5.1)
stringr * 1.6.0 2025-11-04 [1] CRAN (R 4.5.1)
SummarizedExperiment 1.40.0 2025-10-29 [1] Bioconductor 3.22 (R 4.5.1)
survival 3.8-3 2024-12-17 [2] CRAN (R 4.5.1)
tensor 1.5.1 2025-06-17 [1] CRAN (R 4.5.1)
tibble 3.3.0 2025-06-08 [1] CRAN (R 4.5.1)
tidyr 1.3.1 2024-01-24 [1] CRAN (R 4.5.1)
tidyselect 1.2.1 2024-03-11 [1] CRAN (R 4.5.1)
UCSC.utils 1.6.0 2025-10-29 [1] Bioconductor 3.22 (R 4.5.1)
usethis 3.2.1 2025-09-06 [1] CRAN (R 4.5.1)
uwot 0.2.4 2025-11-10 [1] CRAN (R 4.5.1)
VariantAnnotation 1.56.0 2025-10-29 [1] Bioconductor 3.22 (R 4.5.1)
vctrs 0.6.5 2023-12-01 [1] CRAN (R 4.5.1)
vipor 0.4.7 2023-12-18 [1] CRAN (R 4.5.1)
viridisLite 0.4.2 2023-05-02 [1] CRAN (R 4.5.1)
withr 3.0.2 2024-10-28 [1] CRAN (R 4.5.1)
xfun 0.54 2025-10-30 [1] CRAN (R 4.5.1)
XML 3.99-0.20 2025-11-08 [1] CRAN (R 4.5.1)
xtable 1.8-4 2019-04-21 [1] CRAN (R 4.5.1)
XVector 0.50.0 2025-10-29 [1] Bioconductor 3.22 (R 4.5.1)
yaml 2.3.10 2024-07-26 [1] CRAN (R 4.5.1)
zoo 1.8-14 2025-04-10 [1] CRAN (R 4.5.1)
[1] /home/rsg/R/x86_64-pc-linux-gnu-library/4.5
[2] /opt/R/4.5.1/lib/R/library
* ── Packages attached to the search path.
──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────