Analyzing adult mouse brain scATAC-seq
Compiled: December 07, 2022
This vignette echoes the commands run in the introductory Signac vignette on human PBMC. We provide the same analysis in a different system to demonstrate performance and applicability to other tissue types, and to provide an example from another species.
First load in Signac, Seurat, and some other packages we will be using for analyzing mouse data.
library(Signac)
library(Seurat)
library(GenomeInfoDb)
library(EnsDb.Mmusculus.v79)
library(ggplot2)
library(patchwork)
set.seed(1234)Pre-processing workflow
counts <- Read10X_h5("./vignette_data/atac_v1_adult_brain_fresh_5k_filtered_peak_bc_matrix.h5")
metadata <- read.csv(
file = "./vignette_data/atac_v1_adult_brain_fresh_5k_singlecell.csv",
header = TRUE,
row.names = 1
)
brain_assay <- CreateChromatinAssay(
counts = counts,
sep = c(":", "-"),
genome = "mm10",
fragments = './vignette_data/atac_v1_adult_brain_fresh_5k_fragments.tsv.gz',
min.cells = 1
)## Computing hashbrain <- CreateSeuratObject(
counts = brain_assay,
assay = 'peaks',
project = 'ATAC',
meta.data = metadata
)## Warning in CreateSeuratObject.Assay(counts = brain_assay, assay = "peaks", :
## Some cells in meta.data not present in provided counts matrix.## Warning: Keys should be one or more alphanumeric characters followed by an
## underscore, setting key from peaks to peaks_We can also add gene annotations to the brain object for
the mouse genome. This will allow downstream functions to pull the gene
annotation information directly from the object.
# extract gene annotations from EnsDb
annotations <- GetGRangesFromEnsDb(ensdb = EnsDb.Mmusculus.v79)
# change to UCSC style since the data was mapped to hg19
seqlevelsStyle(annotations) <- 'UCSC'
# add the gene information to the object
Annotation(brain) <- annotationsComputing QC Metrics
Next we compute some useful per-cell QC metrics.
brain <- NucleosomeSignal(object = brain)We can 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 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 > 4, 'NS > 4', 'NS < 4')
FragmentHistogram(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.
brain <- TSSEnrichment(brain, fast = FALSE)brain$high.tss <- ifelse(brain$TSS.enrichment > 2, 'High', 'Low')
TSSPlot(brain, group.by = 'high.tss') + NoLegend()
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', 'peak_region_fragments',
'TSS.enrichment', 'blacklist_ratio', 'nucleosome_signal'),
pt.size = 0.1,
ncol = 5
)
We remove cells that are outliers for these QC metrics.
brain <- subset(
x = brain,
subset = peak_region_fragments > 3000 &
peak_region_fragments < 100000 &
pct_reads_in_peaks > 40 &
blacklist_ratio < 0.025 &
nucleosome_signal < 4 &
TSS.enrichment > 2
)
brain## An object of class Seurat
## 276523 features across 3517 samples within 1 assay
## Active assay: peaks (276523 features, 0 variable features)Normalization and linear dimensional reduction
brain <- RunTFIDF(brain)
brain <- FindTopFeatures(brain, min.cutoff = 'q0')
brain <- RunSVD(object = brain)The first LSI component often captures sequencing depth (technical
variation) rather than biological variation. If this is the case, the
component should be removed from downstream analysis. We can assess the
correlation between each LSI component and sequencing depth using the
DepthCor() function:
DepthCor(brain)
Here we see there is a very strong correlation between the first LSI component and the total number of counts for the cell, so we will perform downstream steps without this component.
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 = 2:30
)
brain <- FindNeighbors(
object = brain,
reduction = 'lsi',
dims = 2:30
)
brain <- FindClusters(
object = brain,
algorithm = 3,
resolution = 1.2,
verbose = FALSE
)
DimPlot(object = brain, label = TRUE) + NoLegend()
Create a gene activity matrix
# compute gene activities
gene.activities <- GeneActivity(brain)
# add the gene activity matrix to the Seurat object as a new assay
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("./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 = 2:30
)
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')
plot1 + plot2
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 differentially accessible 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.05,
test.use = 'LR',
latent.vars = 'peak_region_fragments'
)
head(da_peaks)## p_val avg_log2FC pct.1 pct.2 p_val_adj
## chr4-86523460-86525699 7.562228e-68 0.4229695 0.416 0.036 2.091130e-62
## chr4-101303783-101306999 5.082823e-61 0.4036443 0.479 0.074 1.405518e-55
## chr15-87603138-87609228 9.036515e-53 0.3370319 0.569 0.156 2.498804e-47
## chr3-137056226-137059489 1.404579e-51 0.2967401 0.531 0.112 3.883983e-46
## chr19-28380244-28385110 1.035894e-50 0.3210634 0.498 0.104 2.864484e-45
## chr13-69329769-69332192 1.097145e-50 -0.3259721 0.150 0.454 3.033858e-45plot1 <- VlnPlot(
object = brain,
features = rownames(da_peaks)[1],
pt.size = 0.1,
idents = c("L4","L5 IT","L2/3 IT")
)
plot2 <- FeaturePlot(
object = brain,
features = rownames(da_peaks)[1],
pt.size = 0.1,
max.cutoff = 'q95'
)
plot1 | plot2
open_l23 <- rownames(da_peaks[da_peaks$avg_log2FC > 0.25, ])
open_l456 <- rownames(da_peaks[da_peaks$avg_log2FC < -0.25, ])
closest_l23 <- ClosestFeature(brain, open_l23)
closest_l456 <- ClosestFeature(brain, open_l456)
head(closest_l23)## tx_id gene_name gene_id
## ENSMUST00000151481 ENSMUST00000151481 Fam154a ENSMUSG00000028492
## ENSMUST00000131864 ENSMUST00000131864 Gm12796 ENSMUSG00000085721
## ENSMUSE00000647021 ENSMUST00000068088 Fam19a5 ENSMUSG00000054863
## ENSMUST00000070198 ENSMUST00000070198 Ppp3ca ENSMUSG00000028161
## ENSMUST00000162022 ENSMUST00000162022 Glis3 ENSMUSG00000052942
## ENSMUST00000165341 ENSMUST00000165341 Otogl ENSMUSG00000091455
## gene_biotype type closest_region
## ENSMUST00000151481 protein_coding gap chr4-86487920-86538964
## ENSMUST00000131864 lincRNA gap chr4-101292521-101318425
## ENSMUSE00000647021 protein_coding exon chr15-87625230-87625486
## ENSMUST00000070198 protein_coding utr chr3-136935226-136937727
## ENSMUST00000162022 protein_coding gap chr19-28357953-28530872
## ENSMUST00000165341 protein_coding utr chr10-107762223-107762309
## query_region distance
## ENSMUST00000151481 chr4-86523460-86525699 0
## ENSMUST00000131864 chr4-101303783-101306999 0
## ENSMUSE00000647021 chr15-87603138-87609228 16001
## ENSMUST00000070198 chr3-137056226-137059489 118498
## ENSMUST00000162022 chr19-28380244-28385110 0
## ENSMUST00000165341 chr10-107751627-107753416 8806head(closest_l456)## tx_id gene_name gene_id
## ENSMUST00000044081 ENSMUST00000044081 Papd7 ENSMUSG00000034575
## ENSMUSE00000368324 ENSMUST00000027730 Myog ENSMUSG00000026459
## ENSMUST00000189866 ENSMUST00000189866 Iqcf3 ENSMUSG00000023577
## ENSMUST00000076817 ENSMUST00000076817 Utrn ENSMUSG00000019820
## ENSMUST00000037205 ENSMUST00000037205 Mcee ENSMUSG00000033429
## ENSMUST00000097303 ENSMUST00000097303 Arrdc5 ENSMUSG00000073380
## gene_biotype type closest_region
## ENSMUST00000044081 protein_coding utr chr13-69497959-69499915
## ENSMUSE00000368324 protein_coding exon chr1-134289989-134290526
## ENSMUST00000189866 protein_coding gap chr9-106527019-106544398
## ENSMUST00000076817 protein_coding cds chr10-12684410-12684571
## ENSMUST00000037205 protein_coding utr chr7-64411988-64412121
## ENSMUST00000097303 protein_coding cds chr17-56294146-56294664
## query_region distance
## ENSMUST00000044081 chr13-69329769-69332192 165766
## ENSMUSE00000368324 chr1-134272695-134276899 13089
## ENSMUST00000189866 chr9-106534311-106536341 0
## ENSMUST00000076817 chr10-12683407-12684772 0
## ENSMUST00000037205 chr7-64450369-64453260 38247
## ENSMUST00000097303 chr17-56293906-56295053 0Plotting genomic regions
We can also create coverage plots grouped by cluster, cell type, or
any other metadata stored in the object for any genomic region using the
CoveragePlot() function. These represent pseudo-bulk
accessibility tracks, where signal from all cells within a group have
been averaged together to visualize the DNA accessibility in a
region.
# 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")
CoveragePlot(
object = brain,
region = c("Neurod6", "Gad2"),
extend.upstream = 1000,
extend.downstream = 1000,
ncol = 1
)
Session Info
sessionInfo()## R version 4.2.1 (2022-06-23 ucrt)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 19045)
##
## Matrix products: default
##
## locale:
## [1] LC_COLLATE=English_United States.utf8
## [2] LC_CTYPE=English_United States.utf8
## [3] LC_MONETARY=English_United States.utf8
## [4] LC_NUMERIC=C
## [5] LC_TIME=English_United States.utf8
##
## attached base packages:
## [1] stats4 stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] patchwork_1.1.1 ggplot2_3.3.6
## [3] EnsDb.Mmusculus.v79_2.99.0 ensembldb_2.20.2
## [5] AnnotationFilter_1.20.0 GenomicFeatures_1.48.3
## [7] AnnotationDbi_1.58.0 Biobase_2.56.0
## [9] GenomicRanges_1.48.0 GenomeInfoDb_1.32.2
## [11] IRanges_2.30.0 S4Vectors_0.34.0
## [13] BiocGenerics_0.42.0 sp_1.5-0
## [15] SeuratObject_4.1.0 Seurat_4.1.1
## [17] Signac_1.7.0
##
## loaded via a namespace (and not attached):
## [1] utf8_1.2.2 reticulate_1.25
## [3] tidyselect_1.1.2 RSQLite_2.2.14
## [5] htmlwidgets_1.5.4 grid_4.2.1
## [7] BiocParallel_1.30.3 Rtsne_0.16
## [9] munsell_0.5.0 codetools_0.2-18
## [11] ica_1.0-3 interp_1.1-3
## [13] future_1.27.0 miniUI_0.1.1.1
## [15] withr_2.5.0 spatstat.random_2.2-0
## [17] colorspace_2.0-3 progressr_0.10.1
## [19] filelock_1.0.2 highr_0.9
## [21] knitr_1.39 rstudioapi_0.13
## [23] ROCR_1.0-11 tensor_1.5
## [25] listenv_0.8.0 labeling_0.4.2
## [27] MatrixGenerics_1.8.1 GenomeInfoDbData_1.2.8
## [29] polyclip_1.10-0 farver_2.1.1
## [31] bit64_4.0.5 parallelly_1.32.1
## [33] vctrs_0.4.1 generics_0.1.3
## [35] xfun_0.31 biovizBase_1.44.0
## [37] BiocFileCache_2.4.0 R6_2.5.1
## [39] hdf5r_1.3.5 bitops_1.0-7
## [41] spatstat.utils_2.3-1 cachem_1.0.6
## [43] DelayedArray_0.22.0 assertthat_0.2.1
## [45] promises_1.2.0.1 BiocIO_1.6.0
## [47] scales_1.2.0 nnet_7.3-17
## [49] rgeos_0.5-9 gtable_0.3.0
## [51] globals_0.15.1 goftest_1.2-3
## [53] rlang_1.0.3 RcppRoll_0.3.0
## [55] splines_4.2.1 rtracklayer_1.56.1
## [57] lazyeval_0.2.2 dichromat_2.0-0.1
## [59] checkmate_2.1.0 spatstat.geom_2.4-0
## [61] yaml_2.3.5 reshape2_1.4.4
## [63] abind_1.4-5 backports_1.4.1
## [65] httpuv_1.6.5 Hmisc_4.7-0
## [67] tools_4.2.1 ellipsis_0.3.2
## [69] spatstat.core_2.4-4 jquerylib_0.1.4
## [71] RColorBrewer_1.1-3 ggridges_0.5.3
## [73] Rcpp_1.0.8.3 plyr_1.8.7
## [75] base64enc_0.1-3 progress_1.2.2
## [77] zlibbioc_1.42.0 purrr_0.3.4
## [79] RCurl_1.98-1.7 prettyunits_1.1.1
## [81] rpart_4.1.16 deldir_1.0-6
## [83] pbapply_1.5-0 cowplot_1.1.1
## [85] zoo_1.8-10 SummarizedExperiment_1.26.1
## [87] ggrepel_0.9.1 cluster_2.1.3
## [89] magrittr_2.0.3 RSpectra_0.16-1
## [91] data.table_1.14.2 scattermore_0.8
## [93] lmtest_0.9-40 RANN_2.6.1
## [95] ProtGenerics_1.28.0 fitdistrplus_1.1-8
## [97] matrixStats_0.62.0 hms_1.1.1
## [99] mime_0.12 evaluate_0.16
## [101] xtable_1.8-4 XML_3.99-0.10
## [103] jpeg_0.1-9 gridExtra_2.3
## [105] compiler_4.2.1 biomaRt_2.52.0
## [107] tibble_3.1.7 KernSmooth_2.23-20
## [109] crayon_1.5.1 htmltools_0.5.2
## [111] mgcv_1.8-40 later_1.3.0
## [113] Formula_1.2-4 tidyr_1.2.0
## [115] DBI_1.1.3 dbplyr_2.2.1
## [117] MASS_7.3-57 rappdirs_0.3.3
## [119] Matrix_1.4-1 cli_3.3.0
## [121] parallel_4.2.1 igraph_1.3.2
## [123] pkgconfig_2.0.3 GenomicAlignments_1.32.0
## [125] foreign_0.8-82 plotly_4.10.0
## [127] spatstat.sparse_2.1-1 xml2_1.3.3
## [129] bslib_0.4.0 XVector_0.36.0
## [131] VariantAnnotation_1.42.1 stringr_1.4.0
## [133] digest_0.6.29 sctransform_0.3.3
## [135] RcppAnnoy_0.0.19 spatstat.data_2.2-0
## [137] Biostrings_2.64.0 rmarkdown_2.14
## [139] leiden_0.4.2 fastmatch_1.1-3
## [141] htmlTable_2.4.1 uwot_0.1.11
## [143] restfulr_0.0.15 curl_4.3.2
## [145] shiny_1.7.2 Rsamtools_2.12.0
## [147] rjson_0.2.21 lifecycle_1.0.1
## [149] nlme_3.1-157 jsonlite_1.8.0
## [151] BSgenome_1.64.0 viridisLite_0.4.0
## [153] fansi_1.0.3 pillar_1.8.0
## [155] lattice_0.20-45 KEGGREST_1.36.3
## [157] fastmap_1.1.0 httr_1.4.3
## [159] survival_3.3-1 glue_1.6.2
## [161] png_0.1-7 bit_4.0.4
## [163] stringi_1.7.6 sass_0.4.2
## [165] blob_1.2.3 latticeExtra_0.6-30
## [167] memoise_2.0.1 dplyr_1.0.9
## [169] irlba_2.3.5 future.apply_1.9.0