ATAC-seq analysis in R#

In this tutorial, we go over how to use scvi-tools functionality in R for analyzing ATAC-seq data. We will closely follow the PBMC tutorial from Signac, using scvi-tools when appropriate. In particular, we will

  1. Use PeakVI for dimensionality reduction and differential accessiblity for the ATAC-seq data

  2. Use scVI to integrate the unpaired ATAC-seq dataset with a match scRNA-seq dataset of PBMCs

This tutorial requires Reticulate. Please check out our installation guide for instructions on installing Reticulate and scvi-tools.

Loading and processing data with Signac#



We follow the original tutorial to create the Seurat object with ATAC data.

counts <- Read10X_h5(filename = "atac_v1_pbmc_10k_filtered_peak_bc_matrix.h5")
metadata <- read.csv(
  file = "atac_v1_pbmc_10k_singlecell.csv",
  header = TRUE,
  row.names = 1

chrom_assay <- CreateChromatinAssay(
  counts = counts,
  sep = c(":", "-"),
  genome = 'hg19',
  fragments = 'atac_v1_pbmc_10k_fragments.tsv.gz',
  min.cells = 10,
  min.features = 200

pbmc <- CreateSeuratObject(
  counts = chrom_assay,
  assay = "peaks", = metadata
An object of class Seurat 
87561 features across 8728 samples within 1 assay 
Active assay: peaks (87561 features, 0 variable features)
ChromatinAssay data with 87561 features for 8728 cells
Variable features: 0 
Genome: hg19 
Annotation present: FALSE 
Motifs present: FALSE 
Fragment files: 1 

We add gene annotation information to facilitate downstream functionality.

# extract gene annotations from EnsDb
annotations <- GetGRangesFromEnsDb(ensdb = EnsDb.Hsapiens.v75)

# change to UCSC style since the data was mapped to hg19
seqlevelsStyle(annotations) <- 'UCSC'
genome(annotations) <- "hg19"

# add the gene information to the object
Annotation(pbmc) <- annotations

Computing QC metrics#

We compute the same QC metrics as the original tutorial. We leave it to the reader to follow the excellent Signac tutorial for understanding what these quantities represent.

# compute nucleosome signal score per cell
pbmc <- NucleosomeSignal(object = pbmc)

# compute TSS enrichment score per cell
pbmc <- TSSEnrichment(object = pbmc, fast = FALSE)

# add blacklist ratio and fraction of reads in peaks
pbmc$pct_reads_in_peaks <- pbmc$peak_region_fragments / pbmc$passed_filters * 100
pbmc$blacklist_ratio <- pbmc$blacklist_region_fragments / pbmc$peak_region_fragments
options(repr.plot.width=14, repr.plot.height=4)
  object = pbmc,
  features = c('pct_reads_in_peaks', 'peak_region_fragments',
               'TSS.enrichment', 'blacklist_ratio', 'nucleosome_signal'),
  pt.size = 0.1,
  ncol = 5
pbmc <- subset(
  x = pbmc,
  subset = peak_region_fragments > 3000 &
    peak_region_fragments < 20000 &
    pct_reads_in_peaks > 15 &
    blacklist_ratio < 0.05 &
    nucleosome_signal < 4 &
    TSS.enrichment > 2
An object of class Seurat 
87561 features across 7060 samples within 1 assay 
Active assay: peaks (87561 features, 0 variable features)

Dimensionality reduction (PeakVI)#

Creating an AnnData object#

We follow the standard workflow for converting between Seurat and AnnData.

use_python("/home/adam/.pyenv/versions/3.9.7/bin/python", required = TRUE)
sc <- import("scanpy", convert = FALSE)
scvi <- import("scvi", convert = FALSE)
adata <- convertFormat(pbmc, from="seurat", to="anndata", main_layer="counts", assay="peaks", drop_single_values=FALSE)
print(adata) # Note generally in Python, dataset conventions are obs x var
AnnData object with n_obs × n_vars = 7060 × 87561
    obs: 'orig.ident', 'nCount_peaks', 'nFeature_peaks', 'total', 'duplicate', 'chimeric', 'unmapped', 'lowmapq', 'mitochondrial', 'passed_filters', 'cell_id', 'is__cell_barcode', 'TSS_fragments', 'DNase_sensitive_region_fragments', 'enhancer_region_fragments', 'promoter_region_fragments', 'on_target_fragments', 'blacklist_region_fragments', 'peak_region_fragments', 'nucleosome_signal', 'nucleosome_percentile', 'TSS.enrichment', 'TSS.percentile', 'pct_reads_in_peaks', 'blacklist_ratio'
    var: 'count', 'percentile'

Run the standard PeakVI workflow#

pvi <- scvi$model$PEAKVI(adata)
# get the latent represenation
latent = pvi$get_latent_representation()

# put it back in our original Seurat object
latent <- as.matrix(latent)
rownames(latent) = colnames(pbmc)
ndims <- ncol(latent)
pbmc[["peakvi"]] <- CreateDimReducObject(embeddings = latent, key = "peakvi_", assay = "peaks")
# Find clusters, then run UMAP, and visualize
pbmc <- FindNeighbors(pbmc, reduction = "peakvi", dims=1:ndims)
pbmc <- FindClusters(pbmc, resolution = 1)

pbmc <- RunUMAP(pbmc, reduction = "peakvi", dims=1:ndims)

options(repr.plot.width=7, repr.plot.height=7)
DimPlot(object = pbmc, label = TRUE) + NoLegend()
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 7060
Number of edges: 220834

Running Louvain algorithm...
Maximum modularity in 10 random starts: 0.8544
Number of communities: 18
Elapsed time: 0 seconds

Create a gene activity matrix#


The gene activity is used as an approximation of a gene expression matrix such that unpaired ATAC data can be integrated with RNA data. We recommend using this approach only for this unpaired case. Better results can be acheived if there is partially paired data, in which case MultiVI can be used.

gene.activities <- GeneActivity(pbmc)

# add the gene activity matrix to the Seurat object as a new assay and normalize it
pbmc[['RNA']] <- CreateAssayObject(counts = gene.activities)
pbmc <- NormalizeData(
  object = pbmc,
  assay = 'RNA',
  normalization.method = 'LogNormalize',
  scale.factor = median(pbmc$nCount_RNA)
options(repr.plot.width=12, repr.plot.height=10)

DefaultAssay(pbmc) <- 'RNA'

  object = pbmc,
  features = c('MS4A1', 'CD3D', 'LEF1', 'NKG7', 'TREM1', 'LYZ'),
  pt.size = 0.1,
  max.cutoff = 'q95',
  ncol = 3

Integrating with scRNA-seq data (scANVI)#

We can integrate the gene activity matrix with annotated scRNA-seq data using scANVI.

First we download the Seurat-processed PBMC 10k dataset (as in their tutorial).

system("wget -O pbmc_10k_v3.rds")
pbmc_rna <- readRDS("pbmc_10k_v3.rds")

And we convert it to AnnData using sceasy again. Subsequently, we follow the standard scANVI workflow: pretraining with scVI then running scANVI.

adata_rna <- convertFormat(pbmc_rna, from="seurat", to="anndata", main_layer="counts", assay="RNA", drop_single_values=FALSE)
adata_atac_act <- convertFormat(pbmc, from="seurat", to="anndata", main_layer="counts", assay="RNA", drop_single_values=FALSE)

# provide adata_atac_act unknown cell type labels
adata_atac_act$obs$insert(adata_atac_act$obs$shape[1], "celltype", "Unknown")

adata_both <- adata_rna$concatenate(adata_atac_act)

We concatenated the RNA expression with the activity matrix using AnnData. Now we can see the last column is called “batch” and denotes which dataset each cell originated from.

A data.frame: 6 × 35
rna_AAACCCAAGCGCCCAT-1-010x_RNA220410870.0358126720.43820220.023593471CD4 Memory NaNNaNNaNNaNNaNNaNNaNNaNNaNNANA0
rna_AAACCCACAGAGTTGG-1-010x_RNA588418360.0192270340.10179640.107579880CD14+ MonocytesNaNNaNNaNNaNNaNNaNNaNNaNNaNNANA0
rna_AAACCCACAGGTATGG-1-010x_RNA553022160.0054478650.13928010.078481015NK dim NaNNaNNaNNaNNaNNaNNaNNaNNaNNANA0
rna_AAACCCACATAGTCAC-1-010x_RNA510616150.0142760030.49494950.108303963pre-B cell NaNNaNNaNNaNNaNNaNNaNNaNNaNNANA0
rna_AAACCCACATCCAATG-1-010x_RNA457218000.0538573510.13928010.089895015NK bright NaNNaNNaNNaNNaNNaNNaNNaNNaNNANA0
rna_AAACCCAGTGGCTACC-1-010x_RNA670219650.0566037740.35543280.063264701CD4 Memory NaNNaNNaNNaNNaNNaNNaNNaNNaNNANA0
scvi$model$SCVI$setup_anndata(adata_both, labels_key="celltype", batch_key="batch")
model <- scvi$model$SCVI(adata_both, gene_likelihood="nb", dispersion="gene-batch")
lvae <- scvi$model$SCANVI$from_scvi_model(model, "Unknown", adata=adata_both)
lvae$train(max_epochs = as.integer(100), n_samples_per_label = as.integer(100))

Here we only use the prediction functionality of scANVI, but we also could have viewed an integrated representation of the ATAC and RNA using UMAP.

adata_both$obs$insert(adata_both$obs$shape[1], "predicted.labels", lvae$predict())
df <- py_to_r(adata_both$obs)
df <- subset(df, batch == 1)[, c("predicted.labels")]
pbmc <- AddMetaData(object = pbmc, metadata = df,"predicted.labels")


These labels should only serve as a starting point. Further inspection should always be performed. We leave this to the user, but will continue with these labels as a demonstration.

options(repr.plot.width=12, repr.plot.height=5)

plot1 <- DimPlot(
  object = pbmc_rna, = 'celltype',
  label = TRUE,
  repel = TRUE) + NoLegend() + ggtitle('scRNA-seq')

plot2 <- DimPlot(
  object = pbmc, = 'predicted.labels',
  label = TRUE,
  repel = TRUE) + ggtitle('scATAC-seq')

plot1 + plot2

Finding differentially accessible peaks between clusters#

As PeakVI learns uncertainty around the observed data, it can be leveraged for differential accessibility analysis. First, let’s store the seurat cluster information back inside the AnnData.

adata$obs$insert(adata$obs$shape[1], "predicted_ct", pbmc[["predicted.labels"]][,1])

Using our trained PEAKVI model, we call the differential_accessibility() (DA) method We pass predicted_ct to the groupby argument and compare between naive CD4s and CD14 monocytes.

The output of DA is a DataFrame with the bayes factors. Bayes factors > 3 have high probability of being differentially expressed. You can also set fdr_target, which will return the differentially expressed genes based on the posteior expected FDR.

DA <- pvi$differential_accessibility(adata, groupby="predicted_ct", group1 = "CD4 Naive", group2 = "CD14+ Monocytes")
DA <- py_to_r(DA)
A data.frame: 6 × 9
chr1-569174-5696390.5342FALSE 0.1370139-0.06247292-0.0263075850.081191560.0187186390.0355411950.009233610
chr1-713460-7148230.9208 TRUE 2.4532664-0.27283466-0.1640375490.735343690.4625090360.3893376410.225300092
chr1-752422-7530380.9722 TRUE 3.5545252 0.18390924 0.0934914230.018991580.2029008120.0048465270.098337950
chr1-762106-7633590.8202FALSE 1.5177030-0.15394741-0.0966605360.514659700.3607122900.2665589660.169898430
chr1-779589-7802710.7834FALSE 1.2855911-0.09098869-0.0493878820.096615920.0056272280.0516962840.002308403
# sort by proba_da and effect_size
DA <- DA[order(-DA[, 1], -DA[, 4]), ]
A data.frame: 6 × 9
chr17-1510080-15116161.0000TRUE18.420681 0.8669005 0.44714390.042670910.909571410.0274636510.474607572
chr14-50504661-505077770.9998TRUE 8.516943 0.8495569 0.54364860.140019950.989576880.0856219710.629270545
chr19-36397851-364010400.9998TRUE 8.516943 0.8030195 0.39566840.072686460.875705960.0355411950.431209603
chr7-6674416-66753060.9998TRUE 8.516943 0.4002175 0.18144270.016079800.416297260.0096930530.191135734
chr1-23945628-239473090.9998TRUE 8.516943-0.7074655-0.35933360.722085060.014619570.3667205170.007386888
options(repr.plot.width=10, repr.plot.height=5)

DefaultAssay(pbmc) <- 'peaks'
Idents(pbmc) <- pbmc[["predicted.labels"]]

plot1 <- VlnPlot(
  object = pbmc,
  features = rownames(DA)[1],
  pt.size = 0.1,
  idents = c("CD4 Naive","CD14+ Monocytes")
plot2 <- FeaturePlot(
  object = pbmc,
  features = rownames(DA)[1],
  pt.size = 0.1

plot1 | plot2
sI <- sessionInfo()
sI$loadedOnly <- NULL
print(sI, locale=FALSE)
R version 4.1.0 (2021-05-18)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 20.04.2 LTS

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/blas/
LAPACK: /usr/lib/x86_64-linux-gnu/lapack/

attached base packages:
[1] stats4    stats     graphics  grDevices utils     datasets  methods  
[8] base     

other attached packages:
 [1] sceasy_0.0.6              reticulate_1.20          
 [3] patchwork_1.1.1           ggplot2_3.3.5            
 [5] EnsDb.Hsapiens.v75_2.99.0 ensembldb_2.17.4         
 [7] AnnotationFilter_1.17.1   GenomicFeatures_1.45.2   
 [9] AnnotationDbi_1.55.1      Biobase_2.53.0           
[11] GenomicRanges_1.45.0      GenomeInfoDb_1.29.8      
[13] IRanges_2.27.2            S4Vectors_0.31.3         
[15] BiocGenerics_0.39.2       SeuratObject_4.0.2       
[17] Seurat_4.0.4              Signac_1.3.0