For this tutorial, we will be analyzing a single-cell ATAC-seq dataset of adult mouse brain cells provided by 10x Genomics. The following files are used in this vignette, all available through the 10x Genomics website:

This vignette echoes the commands run in the introductory Signac vignette on human PBMC. We provide a secondary analysis in a second system to demonstrate performance on a second system, and to provide a non-human example.

First load in Signac, Seurat, and some other packages we will be using for analyzing human data.

library(Signac)
library(Seurat)
library(GenomeInfoDb)
library(EnsDb.Mmusculus.v79)
library(ggplot2)
set.seed(1234)

Quickstart

The vignette below presents a detailed workflow for the pre-processing, QC, clustering, visualization, integration, and interpretation of single-cell ATAC-seq data. While there are many steps (some of which are computationally intensive), we note that the core analytical components are similar to single-cell RNA-seq, and can be run quickly in a few lines. For example, the following lines encapsulate our scATAC-seq workflow (but does not include the calculation of advanced QC metrics, or integration with scRNA-seq data).

library(dplyr)

counts <- Read10X_h5("/home/stuartt/github/chrom/vignette_data/atac_v1_adult_brain_fresh_5k_filtered_peak_bc_matrix.h5")
metadata <- read.csv("/home/stuartt/github/chrom/vignette_data/atac_v1_adult_brain_fresh_5k_singlecell.csv", row.names = 1)

brain_fast <- CreateSeuratObject(counts = counts,meta.data = metadata) %>%
  subset(peak_region_fragments > 1000) %>%
  RunTFIDF() %>%
  FindTopFeatures(cutoff = 'q75') %>%
  RunSVD(reduction.name='lsi') %>%
  FindNeighbors(reduction='lsi', dims = 1:30) %>%
  FindClusters(resolution = 0.5,verbose = FALSE) %>%
  RunUMAP(reduction='lsi', dims = 1:30)

DimPlot(brain_fast, label = FALSE)

Below, we describe these steps in detail, beginning with pre-processing.

Pre-processing workflow

counts <- Read10X_h5("/home/stuartt/github/chrom/vignette_data/atac_v1_adult_brain_fresh_5k_filtered_peak_bc_matrix.h5")
metadata <- read.csv(
  file = "/home/stuartt/github/chrom/vignette_data/atac_v1_adult_brain_fresh_5k_singlecell.csv",
  header = TRUE,
  row.names = 1
)

brain <- CreateSeuratObject(
  counts = counts,
  assay = 'peaks',
  project = 'ATAC',
  min.cells = 1,
  meta.data = metadata
)
## Warning in CreateSeuratObject(counts = counts, assay = "peaks", project =
## "ATAC", : Some cells in meta.data not present in provided counts matrix.
fragment.path <- '/home/stuartt/github/chrom/vignette_data/atac_v1_adult_brain_fresh_5k_fragments.tsv.gz'

brain <- SetFragments(
  object = brain,
  file = fragment.path
)

Optional step: Filtering the fragment file

To make downstream steps that use this file faster, we can filter the fragments file to contain only reads from cells that we retain in the analysis. This is optional, but slow, and only needs to be performed once. Running this command writes a new file to disk and indexes the file so it is ready to be used by Signac.

fragment_file_filtered <- "/home/stuartt/github/chrom/vignette_data/atac_v1_adult_brain_fresh_5k_fragments_filtered.tsv"
FilterFragments(
  fragment.path = fragment.path
  cells = colnames(brain),
  output.path = fragment_file_filtered
)
brain <- SetFragments(object = brain, file = paste0(fragment_file_filtered, '.bgz'))

Computing QC Metrics

brain <- NucleosomeSignal(object = brain)
brain$pct_reads_in_peaks <- brain$peak_region_fragments / brain$passed_filters * 100
brain$blacklist_ratio <- brain$blacklist_region_fragments / brain$peak_region_fragments

VlnPlot(
  object = brain,
  features = c('pct_reads_in_peaks', 'blacklist_ratio', 'nucleosome_signal', 'peak_region_fragments'),
  pt.size = 0.1,
  ncol = 4) + NoLegend()

We can also look at the fragment length periodicity for all the cells, and group by cells with high or low nucleosomal signal strength. You can see that cells which are outliers for the mononucleosomal/ nucleosome-free ratio (based off the plots above) have different banding patterns. The remaining cells exhibit a pattern that is typical for a successful ATAC-seq experiment.

brain$nucleosome_group <- ifelse(brain$nucleosome_signal > 10, 'NS > 10', 'NS < 10')
PeriodPlot(object = brain, group.by = 'nucleosome_group', region = 'chr1-1-10000000')

The enrichment of Tn5 integration events at transcriptional start sites (TSSs) can also be an important quality control metric to assess the targeting of Tn5 in ATAC-seq experiments. The ENCODE consortium defined a TSS enrichment score as the number of Tn5 integration site around the TSS normalized to the number of Tn5 integration sites in flanking regions. See the ENCODE documentation for more information about the TSS enrichment score (https://www.encodeproject.org/data-standards/terms/). We can calculate the TSS enrichment score for each cell using the TSSEnrichment function in Signac.

# create granges object with TSS positions
gene.ranges <- genes(EnsDb.Mmusculus.v79)
gene.ranges <- gene.ranges[gene.ranges$gene_biotype == 'protein_coding', ]

tss.ranges <- GRanges(
  seqnames = seqnames(gene.ranges),
  ranges = IRanges(start = start(gene.ranges), width = 2),
  strand = strand(gene.ranges)
)

seqlevelsStyle(tss.ranges) <- 'UCSC'
tss.ranges <- keepStandardChromosomes(tss.ranges, pruning.mode = 'coarse')

# to save time use the first 2000 TSSs
brain <- TSSEnrichment(object = brain, tss.positions = tss.ranges[1:2000])
brain$high.tss <- ifelse(brain$TSS.enrichment > 2, 'High', 'Low')
TSSPlot(brain, group.by = 'high.tss') + ggtitle("TSS enrichment score") + NoLegend()

We remove cells that are outliers for these QC metrics.

brain <- subset(brain, subset = peak_region_fragments > 1000 & peak_region_fragments < 100000 & pct_reads_in_peaks > 15 & blacklist_ratio < 0.025 & nucleosome_signal < 10 & TSS.enrichment > 2)
brain
## An object of class Seurat
## 157203 features across 3795 samples within 1 assay
## Active assay: peaks (157203 features)

Normalization and linear dimensional reduction

brain <- RunTFIDF(brain)
brain <- FindTopFeatures(brain, min.cutoff = 'q0')
brain <- RunSVD(
  object = brain,
  assay = 'peaks',
  reduction.key = 'LSI_',
  reduction.name = 'lsi'
)

Non-linear dimension reduction and clustering

Now that the cells are embedded in a low-dimensional space, we can use methods commonly applied for the analysis of scRNA-seq data to perform graph-based clustering, and non-linear dimension reduction for visualization. The functions RunUMAP, FindNeighbors, and FindClusters all come from the Seurat package.

brain <- RunUMAP(
  object = brain,
  reduction = 'lsi',
  dims = 1:30
)
brain <- FindNeighbors(
  object = brain,
  reduction = 'lsi',
  dims = 1:30
)
brain <- FindClusters(
  object = brain,
  verbose = FALSE
)

DimPlot(object = brain, label = TRUE) + NoLegend()

Create a gene activity matrix

#extract gene coordinates from Ensembl, and ensure name formatting is consistent with  Seurat object 
gene.coords <- genes(EnsDb.Mmusculus.v79, filter = ~ gene_biotype == "protein_coding")
seqlevelsStyle(gene.coords) <- 'UCSC'
genebody.coords <- keepStandardChromosomes(gene.coords, pruning.mode = 'coarse')
genebodyandpromoter.coords <- Extend(x = gene.coords, upstream = 2000, downstream = 0)

# build a gene by cell matrix
gene.activities <- FeatureMatrix(
  fragments = fragment.path,
  features = genebodyandpromoter.coords,
  cells = colnames(brain),
  chunk = 10
)

# convert rownames from chromsomal coordinates into gene names
gene.key <- genebodyandpromoter.coords$gene_name
names(gene.key) <- GRangesToString(grange = genebodyandpromoter.coords)
rownames(gene.activities) <- make.unique(gene.key[rownames(gene.activities)])
gene.activities <- gene.activities[rownames(gene.activities)!="",]

#Add the gene activity matrix to the Seurat object as a new assay, and normalize it
brain[['RNA']] <- CreateAssayObject(counts = gene.activities)
brain <- NormalizeData(
  object = brain,
  assay = 'RNA',
  normalization.method = 'LogNormalize',
  scale.factor = median(brain$nCount_RNA)
)
DefaultAssay(brain) <- 'RNA'
FeaturePlot(
  object = brain,
  features = c('Sst','Pvalb',"Gad2","Neurod6","Rorb","Syt6"),
  pt.size = 0.1,
  max.cutoff = 'q95',
  ncol = 3
)

Integrating with scRNA-seq data

To help interpret the scATAC-seq data, we can classify cells based on an scRNA-seq experiment from the same biological system (the adult mouse brain). We utilize methods for cross-modality integration and label transfer, described here, with a more in-depth tutorial here.

You can download the raw data for this experiment from the Allen Institute website, and view the code used to construct this object on GitHub. Alternatively, you can download the pre-processed Seurat object here.

# Load the pre-processed scRNA-seq data
allen_rna <- readRDS("/home/stuartt/github/chrom/vignette_data/allen_brain.rds")
allen_rna <- FindVariableFeatures(
  object = allen_rna,
  nfeatures = 5000
)

transfer.anchors <- FindTransferAnchors(
  reference = allen_rna,
  query = brain,
  reduction = 'cca',
  dims = 1:40
)

predicted.labels <- TransferData(
  anchorset = transfer.anchors,
  refdata = allen_rna$subclass,
  weight.reduction = brain[['lsi']],
  dims = 1:40
)

brain <- AddMetaData(object = brain, metadata = predicted.labels)
plot1 <- DimPlot(allen_rna, group.by = 'subclass', label = TRUE, repel = TRUE) + NoLegend() + ggtitle('scRNA-seq')
plot2 <- DimPlot(brain, group.by = 'predicted.id', label = TRUE, repel = TRUE) + NoLegend() + ggtitle('scATAC-seq')
CombinePlots(list(plot1,plot2), ncol = 2)

Why did we change default parameters?

We changed default parameters for FindIntegrationAnchors and FindVariableFeatures (including more features and dimensions). You can run the analysis both ways, and observe very similar results. However, when using default parameters we mislabel cluster 11 cells as Vip-interneurons, when they are in fact a Meis2 expressing CGE-derived interneuron population recently described by us and others. The reason is that this subset is exceptionally rare in the scRNA-seq data (0.3%), and so the genes define this subset (for example, Meis2) were too lowly expressed to be selected in the initial set of variable features. We therefore need more genes and dimensions to facilitate cross-modality mapping. Interestingly, this subset is 10-fold more abundant in the scATAC-seq data compared to the scRNA-seq data (see this paper for possible explanations.)

You can see that the RNA-based classifications are entirely consistent with the UMAP visualization, computed only on the ATAC-seq data. We can now easily annotate our scATAC-seq derived clusters (alternately, we could use the RNA classifications themselves). We note three small clusters (13, 20, 21) which represent subdivisions of the scRNA-seq labels. Try transferring the ‘cluster’ label (which shows finer distinctions) from the allen scRNA-seq dataset, to annotate them!

# replace each label with its most likely prediction
for(i in levels(brain)) {
  cells_to_reid <- WhichCells(brain, idents = i)
  newid <- names(sort(table(brain$predicted.id[cells_to_reid]),decreasing=TRUE))[1]
  Idents(brain, cells = cells_to_reid) <- newid
}

Find differentially accessible peaks between clusters

Here, we find DA regions between excitatory neurons in different layers of the cortex.

#switch back to working with peaks instead of gene activities
DefaultAssay(brain) <- 'peaks'

da_peaks <- FindMarkers(
  object = brain,
  ident.1 = c("L2/3 IT"),
  ident.2 = c("L4", "L5 IT", "L6 IT"),
  min.pct = 0.4,
  test.use = 'LR',
  latent.vars = 'peak_region_fragments'
)

head(da_peaks)
##                                 p_val  avg_logFC pct.1 pct.2    p_val_adj
## chr2:118700082-118704897 2.796517e-55  0.3249253 0.586 0.173 4.396208e-50
## chr15:87605281-87607659  7.442849e-54  0.4979019 0.452 0.086 1.170038e-48
## chr13:69329933-69331707  9.335275e-54 -0.4557422 0.138 0.433 1.467533e-48
## chr3:137056475-137058371 2.646040e-52  0.3735190 0.498 0.096 4.159655e-47
## chr19:28380552-28383664  3.047770e-48  0.4239796 0.437 0.081 4.791186e-43
## chr7:29070299-29073172   8.313677e-48  0.3662934 0.549 0.156 1.306935e-42
plot1 <- VlnPlot(
  object = brain,
  features = rownames(da_peaks)[1],
  ncol = 3,
  pt.size = 0.1,
  idents = c("L4","L5 IT","L2/3 IT")
)
plot2 <- FeaturePlot(
  object = brain,
  features = rownames(da_peaks)[1],
  ncol = 3,
  pt.size = 0.1,
  max.cutoff = 'q95'
)
CombinePlots(list(plot1,plot2))

open_l23 <- rownames(da_peaks[da_peaks$avg_logFC > 0.25, ])
open_l456 <- rownames(da_peaks[da_peaks$avg_logFC < -0.25, ])
closest_l23 <- ClosestFeature(regions = open_l23, annotation = EnsDb.Mmusculus.v79, sep = c(':', '-'))
closest_l456 <- ClosestFeature(regions = open_l456, annotation = EnsDb.Mmusculus.v79, sep = c(':', '-'))
head(closest_l23)
##                               gene_id gene_name   gene_biotype
## ENSMUSG00000078137 ENSMUSG00000078137   Ankrd63 protein_coding
## ENSMUSG00000054863 ENSMUSG00000054863   Fam19a5 protein_coding
## ENSMUSG00000028161 ENSMUSG00000028161    Ppp3ca protein_coding
## ENSMUSG00000052942 ENSMUSG00000052942     Glis3 protein_coding
## ENSMUSG00000030592 ENSMUSG00000030592      Ryr1 protein_coding
## ENSMUSG00000024998 ENSMUSG00000024998     Plce1 protein_coding
##                    seq_coord_system  symbol entrezid
## ENSMUSG00000078137       chromosome Ankrd63   383787
## ENSMUSG00000054863       chromosome Fam19a5   106014
## ENSMUSG00000028161       chromosome  Ppp3ca    19055
## ENSMUSG00000052942       chromosome   Glis3   226075
## ENSMUSG00000030592       chromosome    Ryr1    20190
## ENSMUSG00000024998       chromosome   Plce1    74055
##                                      region distance
## ENSMUSG00000078137 chr2-118699103-118703963        0
## ENSMUSG00000054863  chr15-87625230-87759336    17570
## ENSMUSG00000028161 chr3-136670124-136937727   118747
## ENSMUSG00000052942  chr19-28258851-28680077        0
## ENSMUSG00000030592   chr7-29003340-29125151        0
## ENSMUSG00000024998  chr19-38481109-38785030        0
head(closest_l456)
##                               gene_id gene_name   gene_biotype
## ENSMUSG00000034575 ENSMUSG00000034575     Papd7 protein_coding
## ENSMUSG00000056222 ENSMUSG00000056222    Spock1 protein_coding
## ENSMUSG00000027827 ENSMUSG00000027827    Kcnab1 protein_coding
## ENSMUSG00000032517 ENSMUSG00000032517      Mobp protein_coding
## ENSMUSG00000031871 ENSMUSG00000031871      Cdh5 protein_coding
## ENSMUSG00000020866 ENSMUSG00000020866   Cacna1g protein_coding
##                    seq_coord_system  symbol entrezid
## ENSMUSG00000034575       chromosome   Papd7   210106
## ENSMUSG00000056222       chromosome  Spock1    20745
## ENSMUSG00000027827       chromosome  Kcnab1    16497
## ENSMUSG00000032517       chromosome    Mobp    17433
## ENSMUSG00000031871       chromosome    Cdh5    12562
## ENSMUSG00000020866       chromosome Cacna1g    12291
##                                      region distance
## ENSMUSG00000034575  chr13-69497959-69533839   166251
## ENSMUSG00000056222  chr13-57421195-57908314        0
## ENSMUSG00000027827   chr3-65109384-65378223        0
## ENSMUSG00000032517 chr9-120149742-120181484        0
## ENSMUSG00000031871 chr8-104101625-104144502        0
## ENSMUSG00000020866  chr11-94408391-94474198        0

We can also create coverage plots grouped by cluster around any genomic region using the CoveragePlot function. These represent ‘pseudo-bulk’ accessibility tracks, where all cells within a cluster have been averaged together, in order to visualize a more robust chromatin landscape.

# set plotting order
levels(brain) <- c("L2/3 IT","L4","L5 IT","L5 PT","L6 CT", "L6 IT","NP","Sst","Pvalb","Vip","Lamp5","Meis2","Oligo","Astro","Endo","VLMC","Macrophage")

region1 <- rownames(da_peaks)[2]
region2 <- GRangesToString(subset(gene.coords, symbol=="Gad2"))

CoveragePlot(
  object = brain,
  region = c(region1, region2),
  sep = c(":", "-"),
  annotation = EnsDb.Mmusculus.v79,
  extend.upstream = 5000,
  extend.downstream = 5000,
  ncol = 1
)

Session Info

## R version 3.6.1 (2019-07-05)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 18.04.3 LTS
##
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.7.1
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.7.1
##
## 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
##
## attached base packages:
## [1] stats4    parallel  stats     graphics  grDevices utils     datasets
## [8] methods   base
##
## other attached packages:
##  [1] dplyr_0.8.3                ggplot2_3.2.1
##  [3] EnsDb.Mmusculus.v79_2.99.0 ensembldb_2.8.1
##  [5] AnnotationFilter_1.8.0     GenomicFeatures_1.36.4
##  [7] AnnotationDbi_1.46.1       Biobase_2.44.0
##  [9] GenomicRanges_1.36.1       GenomeInfoDb_1.20.0
## [11] IRanges_2.18.3             S4Vectors_0.22.1
## [13] BiocGenerics_0.30.0        Seurat_3.1.1.9010
## [15] Signac_0.1.6
##
## loaded via a namespace (and not attached):
##   [1] reticulate_1.13             R.utils_2.9.0
##   [3] tidyselect_0.2.5            poweRlaw_0.70.2
##   [5] RSQLite_2.1.2               htmlwidgets_1.5.1
##   [7] grid_3.6.1                  BiocParallel_1.18.1
##   [9] Rtsne_0.15                  munsell_0.5.0
##  [11] codetools_0.2-16            ica_1.0-2
##  [13] future_1.15.0               withr_2.1.2
##  [15] colorspace_1.4-1            OrganismDbi_1.26.0
##  [17] knitr_1.25                  rstudioapi_0.10
##  [19] ROCR_1.0-7                  gbRd_0.4-11
##  [21] listenv_0.7.0               labeling_0.3
##  [23] Rdpack_0.11-0               GenomeInfoDbData_1.2.1
##  [25] bit64_0.9-7                 rprojroot_1.3-2
##  [27] vctrs_0.2.0                 xfun_0.10
##  [29] biovizBase_1.32.0           ggseqlogo_0.1
##  [31] R6_2.4.0                    rsvd_1.0.2
##  [33] VGAM_1.1-1                  hdf5r_1.3.0
##  [35] bitops_1.0-6                reshape_0.8.8
##  [37] DelayedArray_0.10.0         assertthat_0.2.1
##  [39] SDMTools_1.1-221.1          scales_1.0.0
##  [41] nnet_7.3-12                 gtable_0.3.0
##  [43] npsurv_0.4-0                globals_0.12.4
##  [45] ggbio_1.32.0                seqLogo_1.50.0
##  [47] rlang_0.4.1                 zeallot_0.1.0
##  [49] splines_3.6.1               rtracklayer_1.44.4
##  [51] lazyeval_0.2.2              acepack_1.4.1
##  [53] dichromat_2.0-0             checkmate_1.9.4
##  [55] BiocManager_1.30.9          yaml_2.2.0
##  [57] reshape2_1.4.3              backports_1.1.5
##  [59] Hmisc_4.3-0                 RBGL_1.60.0
##  [61] tools_3.6.1                 ellipsis_0.3.0
##  [63] gplots_3.0.1.1              RColorBrewer_1.1-2
##  [65] ggridges_0.5.1              Rcpp_1.0.3
##  [67] plyr_1.8.4                  base64enc_0.1-3
##  [69] progress_1.2.2              zlibbioc_1.30.0
##  [71] purrr_0.3.3                 RCurl_1.95-4.12
##  [73] prettyunits_1.0.2           rpart_4.1-15
##  [75] pbapply_1.4-2               cowplot_1.0.0
##  [77] zoo_1.8-6                   SummarizedExperiment_1.14.1
##  [79] ggrepel_0.8.1               cluster_2.1.0
##  [81] fs_1.3.1                    motifmatchr_1.6.0
##  [83] magrittr_1.5                RSpectra_0.15-0
##  [85] data.table_1.12.6           lmtest_0.9-37
##  [87] RANN_2.6.1                  ProtGenerics_1.16.0
##  [89] fitdistrplus_1.0-14         matrixStats_0.55.0
##  [91] hms_0.5.2                   lsei_1.2-0
##  [93] evaluate_0.14               xtable_1.8-4
##  [95] XML_3.98-1.20               gridExtra_2.3
##  [97] compiler_3.6.1              biomaRt_2.40.5
##  [99] tibble_2.1.3                KernSmooth_2.23-16
## [101] crayon_1.3.4                R.oo_1.23.0
## [103] htmltools_0.4.0             Formula_1.2-3
## [105] TFBSTools_1.22.0            tidyr_1.0.0
## [107] RcppParallel_4.4.4          DBI_1.0.0
## [109] MASS_7.3-51.4               Matrix_1.2-17
## [111] readr_1.3.1                 R.methodsS3_1.7.1
## [113] gdata_2.18.0                metap_1.1
## [115] igraph_1.2.4.1              pkgconfig_2.0.3
## [117] pkgdown_1.4.1               GenomicAlignments_1.20.1
## [119] TFMPvalue_0.0.8             foreign_0.8-72
## [121] plotly_4.9.1                annotate_1.62.0
## [123] DirichletMultinomial_1.26.0 XVector_0.24.0
## [125] bibtex_0.4.2                VariantAnnotation_1.30.1
## [127] stringr_1.4.0               digest_0.6.22
## [129] sctransform_0.2.0           RcppAnnoy_0.0.14
## [131] tsne_0.1-3                  graph_1.62.0
## [133] CNEr_1.20.0                 Biostrings_2.52.0
## [135] rmarkdown_1.16              leiden_0.3.1
## [137] htmlTable_1.13.2            uwot_0.1.4
## [139] curl_4.2                    Rsamtools_2.0.3
## [141] gtools_3.8.1                lifecycle_0.1.0
## [143] nlme_3.1-142                jsonlite_1.6
## [145] desc_1.2.0                  viridisLite_0.3.0
## [147] BSgenome_1.52.0             pillar_1.4.2
## [149] lattice_0.20-38             GGally_1.4.0
## [151] KEGGREST_1.24.1             httr_1.4.1
## [153] survival_3.1-7              GO.db_3.8.2
## [155] glue_1.3.1                  png_0.1-7
## [157] bit_1.1-14                  stringi_1.4.3
## [159] blob_1.2.0                  latticeExtra_0.6-28
## [161] caTools_1.17.1.2            memoise_1.1.0
## [163] irlba_2.3.3                 future.apply_1.3.0
## [165] ape_5.3