ATAC-seq and Chromatin Accessibility Models

This document covers models for analyzing single-cell ATAC-seq and chromatin accessibility data in scvi-tools.

PeakVI

Purpose: Analysis and integration of single-cell ATAC-seq data using peak counts.

Key Features:

• Variational autoencoder specifically designed for scATAC-seq peak data
• Learns low-dimensional representations of chromatin accessibility
• Performs batch correction across samples
• Enables differential accessibility testing
• Integrates multiple ATAC-seq datasets

When to Use:

• Analyzing scATAC-seq peak count matrices
• Integrating multiple ATAC-seq experiments
• Batch correction of chromatin accessibility data
• Dimensionality reduction for ATAC-seq
• Differential accessibility analysis between cell types or conditions

Data Requirements:

• Peak count matrix (cells × peaks)
• Binary or count data for peak accessibility
• Batch/sample annotations (optional, for batch correction)

Basic Usage:

import scvi

# Prepare data (peaks should be in adata.X)
# Optional: filter peaks
sc.pp.filter_genes(adata, min_cells=3)

# Setup data
scvi.model.PEAKVI.setup_anndata(
    adata,
    batch_key="batch"
)

# Train model
model = scvi.model.PEAKVI(adata)
model.train()

# Get latent representation (batch-corrected)
latent = model.get_latent_representation()
adata.obsm["X_PeakVI"] = latent

# Differential accessibility
da_results = model.differential_accessibility(
    groupby="cell_type",
    group1="TypeA",
    group2="TypeB"
)

Key Parameters:

• n_latent: Dimensionality of latent space (default: 10)
• n_hidden: Number of nodes per hidden layer (default: 128)
• n_layers: Number of hidden layers (default: 1)
• region_factors: Whether to learn region-specific factors (default: True)
• latent_distribution: Distribution for latent space ("normal" or "ln")

Outputs:

• get_latent_representation(): Low-dimensional embeddings for cells
• get_accessibility_estimates(): Normalized accessibility values
• differential_accessibility(): Statistical testing for differential peaks
• get_region_factors(): Peak-specific scaling factors

Best Practices:

1. Filter out low-quality peaks (present in very few cells)
2. Include batch information if integrating multiple samples
3. Use latent representations for clustering and UMAP visualization
4. Consider using region_factors=True for datasets with high technical variation
5. Store latent embeddings in adata.obsm for downstream analysis with scanpy

PoissonVI

Purpose: Quantitative analysis of scATAC-seq fragment counts (more detailed than peak counts).

Key Features:

• Models fragment counts directly (not just peak presence/absence)
• Poisson distribution for count data
• Captures quantitative differences in accessibility
• Enables fine-grained analysis of chromatin state

When to Use:

• Analyzing fragment-level ATAC-seq data
• Need quantitative accessibility measurements
• Higher resolution analysis than binary peak calls
• Investigating gradual changes in chromatin accessibility

Data Requirements:

• Fragment count matrix (cells × genomic regions)
• Count data (not binary)

Basic Usage (PoissonVI lives in scvi.external):

scvi.external.POISSONVI.setup_anndata(
    adata,
    batch_key="batch"
)

model = scvi.external.POISSONVI(adata)
model.train()

# Get results
latent = model.get_latent_representation()
accessibility = model.get_normalized_accessibility()

Key Differences from PeakVI:

• PeakVI: Best for standard peak count matrices, faster
• PoissonVI: Best for quantitative fragment counts, more detailed

When to Choose PoissonVI over PeakVI:

• Working with fragment counts rather than called peaks
• Need to capture quantitative differences
• Have high-quality, high-coverage data
• Interested in subtle accessibility changes

scBasset

Purpose: Deep learning approach to scATAC-seq analysis with interpretability and motif analysis.

Key Features:

• Convolutional neural network (CNN) architecture for sequence-based analysis
• Models raw DNA sequences, not just peak counts
• Enables motif discovery and transcription factor (TF) binding prediction
• Provides interpretable feature importance
• Performs batch correction

When to Use:

• Want to incorporate DNA sequence information
• Interested in TF motif analysis
• Need interpretable models (which sequences drive accessibility)
• Analyzing regulatory elements and TF binding sites
• Predicting accessibility from sequence alone

Data Requirements:

• Peak sequences (extracted from genome)
• Peak accessibility matrix
• Genome reference (for sequence extraction)

Basic Usage (scBasset lives in scvi.external):

# scBasset needs per-peak DNA sequences. Add them to the AnnData first;
# this downloads the genome (once) and stores one-hot codes in adata.varm.
scvi.data.add_dna_sequence(
    adata,
    genome_name="hg38",
    install_genome=True,
)

# Register the per-peak sequence code, then train
scvi.external.SCBASSET.setup_anndata(adata, dna_code_key="dna_code")

model = scvi.external.SCBASSET(adata)
model.train()

# Cell embeddings (low-dimensional latent representation)
latent = model.get_latent_representation()

Key Parameters:

• n_latent: Latent space dimensionality
• conv_layers: Number of convolutional layers
• n_filters: Number of filters per conv layer
• filter_size: Size of convolutional filters

Advanced Features:

• In silico mutagenesis: Predict how sequence changes affect accessibility
• Motif enrichment: Identify enriched TF motifs in accessible regions
• Batch correction: Similar to other scvi-tools models
• Transfer learning: Fine-tune on new datasets

Interpretability Tools:

scBasset learns sequence-aware cell and peak embeddings. Transcription-factor activity is assessed by scoring motif sequences against the trained model rather than calling a single importance function. See the scBasset user guide for the current motif-injection / TF-activity workflow.

# Cell embeddings for clustering / visualization
cell_embedding = model.get_latent_representation()

Model Selection for ATAC-seq

PeakVI

Choose when:

• Standard scATAC-seq analysis workflow
• Have peak count matrices (most common format)
• Need fast, efficient batch correction
• Want straightforward differential accessibility
• Prioritize computational efficiency

Advantages:

• Fast training and inference
• Proven track record for scATAC-seq
• Easy integration with scanpy workflow
• Robust batch correction

PoissonVI

Choose when:

• Have fragment-level count data
• Need quantitative accessibility measures
• Interested in subtle differences
• Have high-coverage, high-quality data

Advantages:

• More detailed quantitative information
• Better for gradient changes
• Appropriate statistical model for counts

scBasset

Choose when:

• Want to incorporate DNA sequence
• Need biological interpretation (motifs, TFs)
• Interested in regulatory mechanisms
• Have computational resources for CNN training
• Want predictive power for new sequences

Advantages:

• Sequence-based, biologically interpretable
• Motif and TF analysis built-in
• Predictive modeling capabilities
• In silico perturbation experiments

Workflow Example: Complete ATAC-seq Analysis

import scvi
import scanpy as sc

# 1. Load and preprocess ATAC-seq data
adata = sc.read_h5ad("atac_data.h5ad")

# 2. Filter low-quality peaks
sc.pp.filter_genes(adata, min_cells=10)

# 3. Setup and train PeakVI
scvi.model.PEAKVI.setup_anndata(
    adata,
    batch_key="sample"
)

model = scvi.model.PEAKVI(adata, n_latent=20)
model.train(max_epochs=400)

# 4. Extract latent representation
latent = model.get_latent_representation()
adata.obsm["X_PeakVI"] = latent

# 5. Downstream analysis
sc.pp.neighbors(adata, use_rep="X_PeakVI")
sc.tl.umap(adata)
sc.tl.leiden(adata, key_added="clusters")

# 6. Differential accessibility
da_results = model.differential_accessibility(
    groupby="clusters",
    group1="0",
    group2="1"
)

# 7. Save model
model.save("peakvi_model")

Integration with Gene Expression (RNA+ATAC)

For paired multimodal data (RNA+ATAC from same cells), use MultiVI instead:

from mudata import MuData

# MultiVI is configured from a MuData object (setup_anndata was removed in v1.3)
mdata = MuData({"rna": rna_adata, "atac": atac_adata})
scvi.model.MULTIVI.setup_mudata(
    mdata,
    batch_key="sample",
    modalities={"rna_layer": "rna", "atac_layer": "atac"},
)

model = scvi.model.MULTIVI(
    mdata,
    n_genes=rna_adata.n_vars,
    n_regions=atac_adata.n_vars,
)
model.train()

# Get joint latent space
latent = model.get_latent_representation()

See models-multimodal.md for more details on multimodal integration.

Best Practices for ATAC-seq Analysis

1. Quality Control:

   • Filter cells with very low or very high peak counts
   • Remove peaks present in very few cells
   • Filter mitochondrial and sex chromosome peaks if needed

2. Batch Correction:

   • Always include batch_key if integrating multiple samples
   • Consider technical covariates (sequencing depth, TSS enrichment)

3. Feature Selection:

   • Unlike RNA-seq, all peaks are often used
   • Consider filtering very rare peaks for efficiency

4. Latent Dimensions:

   • Start with n_latent=10-30 depending on dataset complexity
   • Larger values for more heterogeneous datasets

5. Downstream Analysis:

   • Use latent representations for clustering and visualization
   • Link peaks to genes for regulatory analysis
   • Perform motif enrichment on cluster-specific peaks

6. Computational Considerations:

   • ATAC-seq matrices are often very large (many peaks)
   • Consider downsampling peaks for initial exploration
   • Use GPU acceleration for large datasets