scATAC-seq in R
This page displays the Quarto source; it does not run the analysis.
templates/scatac/guide.qmd · Download .qmd · GitHub source
---
title: "scATAC-seq with Signac"
format:
html:
toc: true
code-fold: false
execute:
warning: false
message: false
---
> Status: draft. Match the genome build across peaks, fragments, annotations, and sequence resources.
Input: a peak-by-cell count matrix, cell metadata, fragment file, and gene
annotation. The cells run fragment-aware QC, LSI, exploratory peak markers, and
motif enrichment; sample-level peak testing and chromVAR are optional.
## 1. Create the Signac object
The count matrix has peaks in rows and cells in columns. Metadata has one row
per cell and contains `cell_id` and `sample_id`.
```{r}
library(Seurat)
library(Signac)
library(qs2)
library(readr)
counts <- qs2::qs_read("input/peak_counts.qs2")
metadata <- readr::read_tsv("input/cell_metadata.tsv", show_col_types = FALSE)
annotation <- qs2::qs_read("input/gene_annotation.qs2")
fragments_path <- "input/fragments.tsv.gz"
stopifnot("cell_id" %in% names(metadata), "sample_id" %in% names(metadata))
stopifnot(!anyNA(metadata$cell_id), !anyDuplicated(metadata$cell_id))
stopifnot(file.exists(fragments_path), file.exists(paste0(fragments_path, ".tbi")))
stopifnot(!anyDuplicated(rownames(counts)), !anyDuplicated(colnames(counts)))
cell_alignment <- data.frame(
cell_id = c(setdiff(colnames(counts), metadata$cell_id),
setdiff(metadata$cell_id, colnames(counts))),
source = c(rep("count_matrix_only", length(setdiff(colnames(counts), metadata$cell_id))),
rep("metadata_only", length(setdiff(metadata$cell_id, colnames(counts)))))
)
cells <- intersect(colnames(counts), metadata$cell_id)
counts <- counts[, cells, drop = FALSE]
metadata <- metadata[match(cells, metadata$cell_id), , drop = FALSE]
rownames(metadata) <- metadata$cell_id
fragment_object <- Signac::CreateFragmentObject(
path = fragments_path, cells = cells, validate.fragments = TRUE
)
chromatin_assay <- Signac::CreateChromatinAssay(
counts = counts, sep = c(":", "-"), genome = "hg38",
annotation = annotation, fragments = fragment_object,
min.cells = 1, min.features = 1
)
object <- Seurat::CreateSeuratObject(counts = chromatin_assay,
assay = "ATAC", meta.data = metadata)
DefaultAssay(object) <- "ATAC"
print(object)
```
Change `hg38` and the peak separator to match the input. Fragment files must
be indexed and contain the retained cell barcodes.
## 2. Add fragment-aware QC and filter
`metrics.tsv` has `cell_id`, `passed_filters`, `peak_region_fragments`, and
`blacklist_region_fragments`. These example thresholds are not universal.
```{r}
metrics <- readr::read_tsv("input/singlecell_metrics.tsv", show_col_types = FALSE)
stopifnot(!anyDuplicated(metrics$cell_id))
metrics <- metrics[match(colnames(object), metrics$cell_id), , drop = FALSE]
stopifnot(identical(metrics$cell_id, colnames(object)))
metrics$pct_reads_in_peaks <- 100 * metrics$peak_region_fragments / metrics$passed_filters
metrics$blacklist_ratio <- metrics$blacklist_region_fragments / metrics$peak_region_fragments
rownames(metrics) <- metrics$cell_id
object <- Seurat::AddMetaData(object, metrics[, setdiff(names(metrics), "cell_id"), drop = FALSE])
object <- Signac::NucleosomeSignal(object)
object <- Signac::TSSEnrichment(object, fast = TRUE)
qc_metadata <- object[[]]
object$qc_pass <- with(qc_metadata,
is.finite(nCount_ATAC) & is.finite(passed_filters) &
is.finite(pct_reads_in_peaks) & is.finite(blacklist_ratio) &
is.finite(nucleosome_signal) & is.finite(TSS.enrichment) &
nCount_ATAC >= 1000 & passed_filters >= 3000 &
pct_reads_in_peaks >= 25 & blacklist_ratio <= 0.05 &
nucleosome_signal <= 1.5 & TSS.enrichment >= 6
)
print(table(object$sample_id, object$qc_pass, useNA = "ifany"))
VlnPlot(object, features = c("nCount_ATAC", "nFeature_ATAC",
"nucleosome_signal", "TSS.enrichment"),
group.by = "sample_id", ncol = 2, pt.size = 0)
object <- subset(object, subset = qc_pass)
```
Inspect missing values and distributions before filtering. Retain the
per-sample counts so the filter cannot hide a library-specific failure.
## 3. Run TF-IDF, LSI, clustering, and UMAP
```{r}
object <- Signac::RunTFIDF(object)
object <- Signac::FindTopFeatures(object, min.cutoff = "q5")
object <- Signac::RunSVD(object, n = 50)
Signac::DepthCor(object, reduction = "lsi")
# Inspect DepthCor before deciding whether to include LSI1.
object <- Seurat::FindNeighbors(object, reduction = "lsi", dims = 2:30,
graph.name = c("atac_nn", "atac_snn"))
object <- Seurat::FindClusters(object, graph.name = "atac_snn",
resolution = 0.5, random.seed = 1)
object <- Seurat::RunUMAP(object, reduction = "lsi", dims = 2:30,
reduction.name = "umap_atac", seed.use = 1)
DimPlot(object, reduction = "umap_atac", group.by = "seurat_clusters", label = TRUE)
```
LSI is a representation of accessibility. It does not correct batch effects.
## 4. Find differential peaks and motif enrichment
This cell-level cluster comparison is exploratory. `nCount_ATAC` adjusts for
depth, but cells from the same sample are not independent biological
replicates. Use the sample-level pseudobulk section for condition-level claims.
The motif results below inherit that exploratory foreground. The background is
all peaks in the assay; change it only when the study has a justified
background universe.
```{r}
library(JASPAR2024)
library(TFBSTools)
library(RSQLite)
library(BSgenome.Hsapiens.UCSC.hg38)
Seurat::Idents(object) <- "seurat_clusters"
peak_results <- Seurat::FindMarkers(object, ident.1 = "0", ident.2 = "1",
assay = "ATAC", test.use = "LR",
latent.vars = "nCount_ATAC")
foreground <- rownames(peak_results)[peak_results$p_val_adj < 0.05]
background <- rownames(object[["ATAC"]])
stopifnot(length(foreground) > 0)
database <- JASPAR2024::JASPAR2024()
connection <- RSQLite::dbConnect(RSQLite::SQLite(), JASPAR2024::db(database))
motifs <- TFBSTools::getMatrixSet(connection, opts = list(
collection = "CORE", tax_group = "vertebrates", species = 9606,
all_versions = FALSE
))
RSQLite::dbDisconnect(connection)
object <- Signac::AddMotifs(object, genome = BSgenome.Hsapiens.UCSC.hg38,
pfm = motifs, assay = "ATAC")
motif_results <- Signac::FindMotifs(object, features = foreground,
background = background, assay = "ATAC")
head(motif_results[order(motif_results$p.adjust), ], 20)
```
## 5. Test condition differences by sample (optional)
Use one row per biological sample in `sample_metadata.tsv`, with `sample_id`
and `condition` plus any prespecified sample-level covariates. This sums peak
counts within each sample and cell type before testing. The example expects
`control` and `case`; edit the factor levels and contrast for other labels. Add
only estimable covariates to the design; a full-rank design cannot resolve
condition-batch confounding.
```{r}
#| eval: false
library(edgeR)
library(limma)
sample_metadata <- readr::read_tsv("input/sample_metadata.tsv", show_col_types = FALSE)
stopifnot(all(c("sample_id", "condition") %in% names(sample_metadata)))
stopifnot(!anyDuplicated(sample_metadata$sample_id))
stopifnot(!anyNA(sample_metadata$sample_id), !anyNA(sample_metadata$condition))
stopifnot(all(unique(object$sample_id) %in% sample_metadata$sample_id))
sample_covariates <- setdiff(names(sample_metadata), c("sample_id", "condition"))
sample_metadata <- sample_metadata[match(object$sample_id, sample_metadata$sample_id), , drop = FALSE]
rownames(sample_metadata) <- colnames(object)
object <- Seurat::AddMetaData(
object,
sample_metadata[, setdiff(names(sample_metadata), "sample_id"), drop = FALSE]
)
object$cell_type <- as.character(object$seurat_clusters)
print(table(object$sample_id, object$cell_type))
pseudobulk <- Seurat::AggregateExpression(
object,
assays = "ATAC",
return.seurat = TRUE,
group.by = c("sample_id", "condition", "cell_type", sample_covariates)
)
cell_type_to_test <- "0"
pseudobulk <- subset(pseudobulk, subset = cell_type == cell_type_to_test)
pb_counts <- as.matrix(Seurat::GetAssayData(pseudobulk, assay = "ATAC", layer = "counts"))
pb_metadata <- pseudobulk[[]][colnames(pb_counts), , drop = FALSE]
stopifnot(all(as.character(pb_metadata$condition) %in% c("control", "case")))
pb_metadata$condition <- factor(pb_metadata$condition, levels = c("control", "case"))
print(table(pb_metadata$condition))
stopifnot(!anyDuplicated(pb_metadata$sample_id))
stopifnot(all(table(pb_metadata$condition) >= 2))
model_covariates <- intersect("batch", sample_covariates)
# For paired samples, add "donor_id" only when the design is paired:
# model_covariates <- c(model_covariates, "donor_id")
design <- stats::model.matrix(
stats::reformulate(c("condition", model_covariates), intercept = FALSE),
data = pb_metadata
)
stopifnot(qr(design)$rank == ncol(design))
dge <- edgeR::DGEList(counts = pb_counts)
keep <- edgeR::filterByExpr(dge, design = design)
stopifnot(any(keep))
dge <- dge[keep, , keep.lib.sizes = FALSE]
dge <- edgeR::normLibSizes(dge)
fit <- edgeR::glmQLFit(dge, design)
contrast <- limma::makeContrasts(conditioncase - conditioncontrol, levels = design)
test <- edgeR::glmQLFTest(fit, contrast = contrast)
pseudobulk_peak_results <- edgeR::topTags(test, n = nrow(dge))$table
pseudobulk_peak_results$peak <- rownames(pseudobulk_peak_results)
pseudobulk_foreground <- rownames(pseudobulk_peak_results)[pseudobulk_peak_results$FDR < 0.05]
dir.create("output", recursive = TRUE, showWarnings = FALSE)
readr::write_tsv(pseudobulk_peak_results, "output/pseudobulk_differential_peaks.tsv")
if (length(pseudobulk_foreground) > 0) {
population_motif_results <- Signac::FindMotifs(
object, features = pseudobulk_foreground, background = background, assay = "ATAC"
)
readr::write_tsv(population_motif_results, "output/pseudobulk_motif_enrichment.tsv")
}
```
Two samples per condition is only a minimal code check, not evidence of
adequate replication. The disabled block reruns motif enrichment with the
pseudobulk peak foreground.
## 6. Estimate cell-level motif deviations with chromVAR (optional)
Motif enrichment and chromVAR answer different questions. The first tests a
peak set; chromVAR estimates bias-corrected deviations per cell.
```{r}
#| eval: false
library(chromVAR)
peak_counts <- Seurat::GetAssayData(object, assay = "ATAC", layer = "counts")
motif_annotations <- Signac::GetMotifData(object, assay = "ATAC", slot = "data")
deviations <- chromVAR::computeDeviations(
object = peak_counts,
annotations = t(motif_annotations)
)
motif_deviations <- chromVAR::deviationScores(deviations)
dir.create("output", recursive = TRUE, showWarnings = FALSE)
readr::write_tsv(data.frame(motif = rownames(motif_deviations), motif_deviations),
"output/chromvar_deviations.tsv")
```
## 7. Save the object and peak results
```{r}
dir.create("output", recursive = TRUE, showWarnings = FALSE)
qs2::qs_save(object, "output/scatac_analysis.qs2")
readr::write_tsv(data.frame(peak = rownames(peak_results), peak_results),
"output/cell_level_cluster_peak_markers.tsv")
readr::write_tsv(data.frame(motif = rownames(motif_results), motif_results),
"output/motif_enrichment.tsv")
readr::write_tsv(cell_alignment, "output/cell_alignment.tsv")
```