scATAC-seq Analysis with PeakVI

This reference covers single-cell ATAC-seq analysis using PeakVI for dimensionality reduction, batch correction, and differential accessibility.

Overview

PeakVI is a deep generative model for scATAC-seq data that:

• Models binary accessibility (peak open/closed)
• Handles batch effects
• Provides latent representation for clustering
• Enables differential accessibility analysis

Prerequisites

import scvi
import scanpy as sc
import numpy as np
import anndata as ad

print(f"scvi-tools version: {scvi.__version__}")

Step 1: Load and Prepare ATAC Data

From 10x Genomics (Cell Ranger ATAC)

# Peak-cell matrix from fragments
# Usually in filtered_peak_bc_matrix format

adata = sc.read_10x_h5("filtered_peak_bc_matrix.h5")

# Or from mtx format
adata = sc.read_10x_mtx("filtered_peak_bc_matrix/")

# Check structure
print(f"Cells: {adata.n_obs}, Peaks: {adata.n_vars}")
print(f"Sparsity: {1 - adata.X.nnz / (adata.n_obs * adata.n_vars):.2%}")

From ArchR/Signac

# Export from ArchR (in R)
# saveArchRProject(proj, outputDirectory="atac_export", load=FALSE)
# Then read the exported files in Python

# From Signac:
# Export peak matrix and metadata

Step 2: Quality Control

# Calculate QC metrics
sc.pp.calculate_qc_metrics(adata, inplace=True)

# Key metrics for ATAC:
# - n_genes_by_counts: peaks per cell (should rename)
# - total_counts: fragments per cell
adata.obs['n_peaks'] = adata.obs['n_genes_by_counts']
adata.obs['total_fragments'] = adata.obs['total_counts']

# Filter cells
adata = adata[adata.obs['n_peaks'] > 500].copy()
adata = adata[adata.obs['n_peaks'] < 50000].copy()  # Remove potential doublets

# Filter peaks (accessible in at least n cells)
sc.pp.filter_genes(adata, min_cells=10)

print(f"After QC: {adata.shape}")

Binarize Data

# PeakVI works with binary accessibility
# Binarize if not already binary
adata.X = (adata.X > 0).astype(np.float32)

# Verify
print(f"Unique values: {np.unique(adata.X.data)}")

Step 3: Feature Selection

Unlike RNA-seq, peak selection for ATAC is less established. Options:

Option A: Most Accessible Peaks

# Select top peaks by accessibility frequency
peak_accessibility = np.array(adata.X.sum(axis=0)).flatten()
top_peaks = np.argsort(peak_accessibility)[-50000:]  # Top 50k peaks

adata = adata[:, top_peaks].copy()

Option B: Variable Peaks

# Select peaks with high variance
# (Most informative for clustering)
from sklearn.feature_selection import VarianceThreshold

selector = VarianceThreshold(threshold=0.05)
selector.fit(adata.X)
adata = adata[:, selector.get_support()].copy()

Option C: Peaks Near Genes

# Keep peaks within promoter regions or gene bodies
# Requires peak annotation
# gene_peaks = peaks with gene annotation
# adata = adata[:, adata.var['near_gene']].copy()

Step 4: Add Batch Information

# Add batch annotation if multiple samples
adata.obs['batch'] = adata.obs['sample_id']  # Or appropriate column

print(adata.obs['batch'].value_counts())

Step 5: Setup and Train PeakVI

# Setup AnnData
scvi.model.PEAKVI.setup_anndata(
    adata,
    batch_key="batch"  # Optional, omit for single batch
)

# Create model
model = scvi.model.PEAKVI(
    adata,
    n_latent=20,      # Latent dimensions
    n_layers_encoder=2,
    n_layers_decoder=2
)

# Train
model.train(
    max_epochs=200,
    early_stopping=True,
    batch_size=128
)

# Check training
model.history['elbo_train'].plot()

Step 6: Get Latent Representation

# Latent space for downstream analysis
adata.obsm["X_PeakVI"] = model.get_latent_representation()

# Clustering and visualization
sc.pp.neighbors(adata, use_rep="X_PeakVI", n_neighbors=15)
sc.tl.umap(adata)
sc.tl.leiden(adata, resolution=0.5)

# Visualize
sc.pl.umap(adata, color=['leiden', 'batch'], ncols=2)

Step 7: Differential Accessibility

# Differential accessibility between clusters
da_results = model.differential_accessibility(
    groupby='leiden',
    group1='0',
    group2='1'
)

# Filter significant peaks
da_sig = da_results[
    (da_results['is_da_fdr_0.05']) &
    (abs(da_results['lfc_mean']) > 1)
]

print(f"Significant DA peaks: {len(da_sig)}")
print(da_sig.head())

DA Between Conditions

# Compare conditions within cell type
adata_subset = adata[adata.obs['cell_type'] == 'CD4 T cells'].copy()

da_condition = model.differential_accessibility(
    groupby='condition',
    group1='treated',
    group2='control'
)

Step 8: Peak Annotation

# Annotate peaks with nearest genes
# Using pybedtools or similar

# Example peak name format: chr1:1000-2000
# Parse into bed format for annotation

import pandas as pd

def parse_peak_names(peak_names):
    """Parse peak names into bed format."""
    records = []
    for peak in peak_names:
        chrom, coords = peak.split(':')
        start, end = coords.split('-')
        records.append({
            'chrom': chrom,
            'start': int(start),
            'end': int(end),
            'peak': peak
        })
    return pd.DataFrame(records)

peak_bed = parse_peak_names(adata.var_names)

Step 9: Motif Analysis

# Export significant peaks for motif analysis
# Use HOMER, MEME, or chromVAR

# Export peak sequences
sig_peaks = da_sig.index.tolist()
peak_bed_sig = peak_bed[peak_bed['peak'].isin(sig_peaks)]
peak_bed_sig.to_csv("significant_peaks.bed", sep='\t', index=False, header=False)

# Then run HOMER:
# findMotifsGenome.pl significant_peaks.bed hg38 motif_output/ -size 200

Step 10: Gene Activity Scores

# Compute gene activity from peak accessibility
# (Requires peak-gene annotations)

def compute_gene_activity(adata, peak_gene_map):
    """
    Compute gene activity scores from peak accessibility.
    
    Parameters
    ----------
    adata : AnnData
        ATAC data with peaks
    peak_gene_map : dict
        Mapping of peaks to genes
        
    Returns
    -------
    AnnData with gene activity scores
    """
    from scipy.sparse import csr_matrix
    
    genes = list(set(peak_gene_map.values()))
    gene_matrix = np.zeros((adata.n_obs, len(genes)))
    
    for i, gene in enumerate(genes):
        gene_peaks = [p for p, g in peak_gene_map.items() if g == gene]
        if gene_peaks:
            peak_idx = [list(adata.var_names).index(p) for p in gene_peaks if p in adata.var_names]
            if peak_idx:
                gene_matrix[:, i] = np.array(adata.X[:, peak_idx].sum(axis=1)).flatten()
    
    adata_gene = ad.AnnData(
        X=csr_matrix(gene_matrix),
        obs=adata.obs.copy(),
        var=pd.DataFrame(index=genes)
    )
    
    return adata_gene

Complete Pipeline

def analyze_scatac(
    adata,
    batch_key=None,
    n_top_peaks=50000,
    n_latent=20,
    resolution=0.5
):
    """
    Complete scATAC-seq analysis with PeakVI.
    
    Parameters
    ----------
    adata : AnnData
        Raw peak-cell matrix
    batch_key : str, optional
        Batch annotation column
    n_top_peaks : int
        Number of top peaks to use
    n_latent : int
        Latent dimensions
    resolution : float
        Leiden clustering resolution
        
    Returns
    -------
    Tuple of (processed AnnData, trained model)
    """
    import scvi
    import scanpy as sc
    import numpy as np
    
    adata = adata.copy()
    
    # QC
    sc.pp.calculate_qc_metrics(adata, inplace=True)
    adata = adata[adata.obs['n_genes_by_counts'] > 500].copy()
    sc.pp.filter_genes(adata, min_cells=10)
    
    # Binarize
    adata.X = (adata.X > 0).astype(np.float32)
    
    # Select top peaks
    if adata.n_vars > n_top_peaks:
        peak_accessibility = np.array(adata.X.sum(axis=0)).flatten()
        top_peaks = np.argsort(peak_accessibility)[-n_top_peaks:]
        adata = adata[:, top_peaks].copy()
    
    # Setup PeakVI
    scvi.model.PEAKVI.setup_anndata(adata, batch_key=batch_key)
    
    # Train
    model = scvi.model.PEAKVI(adata, n_latent=n_latent)
    model.train(max_epochs=200, early_stopping=True)
    
    # Latent representation
    adata.obsm["X_PeakVI"] = model.get_latent_representation()
    
    # Clustering
    sc.pp.neighbors(adata, use_rep="X_PeakVI")
    sc.tl.umap(adata)
    sc.tl.leiden(adata, resolution=resolution)
    
    return adata, model

# Usage
adata, model = analyze_scatac(
    adata,
    batch_key="sample",
    n_top_peaks=50000
)

# Visualize
sc.pl.umap(adata, color=['leiden', 'sample'])

# Differential accessibility
da_results = model.differential_accessibility(
    groupby='leiden',
    group1='0',
    group2='1'
)

Integration with scRNA-seq

For multiome data or separate RNA/ATAC from same cells:

# See MultiVI for joint RNA+ATAC analysis
# Or use WNN (weighted nearest neighbors) approach

# Transfer labels from RNA to ATAC using shared latent space

Troubleshooting

Key References

• Ashuach et al. (2022) "PeakVI: A deep generative model for single-cell chromatin accessibility analysis"

Advanced Batch Correction with sysVI

This reference covers system-level batch correction using sysVI, designed for integrating data across major technological or study differences.

Overview

sysVI (System Variational Inference) extends scVI for scenarios where:

• Batch effects are very strong (different technologies)
• Standard scVI over-corrects biological signal
• You need to separate "system" effects from biological variation

When to Use sysVI vs scVI

Prerequisites

import scvi
import scanpy as sc
import numpy as np

print(f"scvi-tools version: {scvi.__version__}")

Understanding sysVI Architecture

sysVI separates variation into:

1. Biological variation: Cell type, state, trajectory
2. System variation: Technology, study, lab effects

                    ┌─────────────────┐
Input counts ──────►│    Encoder      │
                    │                 │
System info ───────►│  (conditioned)  │
                    └────────┬────────┘
                             │
                    ┌────────▼────────┐
                    │   Latent z      │
                    │  (biological)   │
                    └────────┬────────┘
                             │
                    ┌────────▼────────┐
System info ───────►│    Decoder      │
                    │  (conditioned)  │
                    └────────┬────────┘
                             │
                    Reconstructed counts

Basic sysVI Workflow

Step 1: Prepare Data

# Load datasets from different systems
adata1 = sc.read_h5ad("10x_data.h5ad")
adata2 = sc.read_h5ad("smartseq_data.h5ad")

# Add system labels
adata1.obs["system"] = "10x"
adata2.obs["system"] = "Smart-seq2"

# Add batch labels (within system)
# e.g., different samples within each technology

# Concatenate
adata = sc.concat([adata1, adata2])

# Store raw counts
adata.layers["counts"] = adata.X.copy()

Step 2: HVG Selection

# Select HVGs considering both batch and system
sc.pp.highly_variable_genes(
    adata,
    n_top_genes=4000,  # More genes for cross-system
    flavor="seurat_v3",
    batch_key="system",  # Consider system for HVG
    layer="counts"
)

# Optionally: ensure overlap between systems
# Check HVGs are expressed in both systems
adata = adata[:, adata.var["highly_variable"]].copy()

Step 3: Setup and Train sysVI

# Setup AnnData
# Note: sysVI may be accessed differently depending on version
# Check scvi-tools documentation for current API

scvi.model.SCVI.setup_anndata(
    adata,
    layer="counts",
    batch_key="sample",           # Within-system batches
    categorical_covariate_keys=["system"]  # System-level covariate
)

# For true sysVI (if available in your version)
# scvi.model.SysVI.setup_anndata(...)

# Create model with system awareness
model = scvi.model.SCVI(
    adata,
    n_latent=30,
    n_layers=2,
    gene_likelihood="nb"
)

# Train
model.train(max_epochs=300)

Step 4: Extract Representations

# Get latent representation
adata.obsm["X_integrated"] = model.get_latent_representation()

# Clustering and visualization
sc.pp.neighbors(adata, use_rep="X_integrated")
sc.tl.umap(adata)
sc.tl.leiden(adata)

# Check integration
sc.pl.umap(adata, color=["system", "leiden", "cell_type"])

Alternative: Harmony + scVI

For cross-system integration, combining methods can work well:

import scanpy.external as sce

# First run PCA
sc.pp.pca(adata)

# Apply Harmony for initial alignment
sce.pp.harmony_integrate(adata, key="system")

# Then train scVI on Harmony-corrected embedding
# Or use Harmony representation directly

Alternative: Using Covariates in scVI

For moderate system effects:

# Include system as categorical covariate
scvi.model.SCVI.setup_anndata(
    adata,
    layer="counts",
    batch_key="sample",
    categorical_covariate_keys=["system", "technology_version"]
)

model = scvi.model.SCVI(adata, n_latent=30)
model.train()

Alternative: Separate Models + Integration

For very different systems:

# Train separate models
scvi.model.SCVI.setup_anndata(adata1, layer="counts", batch_key="sample")
model1 = scvi.model.SCVI(adata1)
model1.train()

scvi.model.SCVI.setup_anndata(adata2, layer="counts", batch_key="sample")
model2 = scvi.model.SCVI(adata2)
model2.train()

# Get latent spaces
adata1.obsm["X_scVI"] = model1.get_latent_representation()
adata2.obsm["X_scVI"] = model2.get_latent_representation()

# Align with CCA or Harmony
# ... additional alignment step

Evaluating Cross-System Integration

Visual Assessment

import matplotlib.pyplot as plt

fig, axes = plt.subplots(1, 3, figsize=(15, 4))

# Color by system
sc.pl.umap(adata, color="system", ax=axes[0], show=False, title="By System")

# Color by cell type
sc.pl.umap(adata, color="cell_type", ax=axes[1], show=False, title="By Cell Type")

# Color by expression of marker
sc.pl.umap(adata, color="CD3D", ax=axes[2], show=False, title="CD3D Expression")

plt.tight_layout()

Quantitative Metrics

# Using scib-metrics
from scib_metrics.benchmark import Benchmarker

bm = Benchmarker(
    adata,
    batch_key="system",
    label_key="cell_type",
    embedding_obsm_keys=["X_integrated"]
)

bm.benchmark()

# Key metrics:
# - Batch mixing (ASW_batch, Graph connectivity)
# - Bio conservation (NMI, ARI, ASW_label)

LISI Scores

# Local Inverse Simpson's Index
from scib_metrics import lisi

# Batch LISI (higher = better mixing)
batch_lisi = lisi.ilisi_graph(
    adata,
    batch_key="system",
    use_rep="X_integrated"
)

# Cell type LISI (lower = better preservation)
ct_lisi = lisi.clisi_graph(
    adata,
    label_key="cell_type", 
    use_rep="X_integrated"
)

print(f"Batch LISI: {batch_lisi.mean():.3f}")
print(f"Cell type LISI: {ct_lisi.mean():.3f}")

Handling Specific Challenges

Different Gene Sets

# Find common genes
common_genes = adata1.var_names.intersection(adata2.var_names)
print(f"Common genes: {len(common_genes)}")

# If too few, use gene mapping
# Or impute missing genes

Different Sequencing Depths

# Add depth as continuous covariate
adata.obs["log_counts"] = np.log1p(adata.obs["total_counts"])

scvi.model.SCVI.setup_anndata(
    adata,
    layer="counts",
    batch_key="sample",
    continuous_covariate_keys=["log_counts"]
)

Unbalanced Cell Types

# Check cell type distribution per system
import pandas as pd

ct_dist = pd.crosstab(adata.obs["system"], adata.obs["cell_type"], normalize="index")
print(ct_dist)

# If very unbalanced, consider:
# 1. Subsample to balance
# 2. Use scANVI with labels to preserve rare types

Complete Pipeline

def integrate_cross_system(
    adatas: dict,
    system_key: str = "system",
    batch_key: str = "batch",
    cell_type_key: str = "cell_type",
    n_top_genes: int = 4000,
    n_latent: int = 30
):
    """
    Integrate datasets from different technological systems.
    
    Parameters
    ----------
    adatas : dict
        Dictionary of {system_name: AnnData}
    system_key : str
        Key for system annotation
    batch_key : str
        Key for within-system batch
    cell_type_key : str
        Key for cell type labels (optional)
    n_top_genes : int
        Number of HVGs
    n_latent : int
        Latent dimensions
        
    Returns
    -------
    Integrated AnnData with model
    """
    import scvi
    import scanpy as sc
    
    # Add system labels and concatenate
    for system_name, adata in adatas.items():
        adata.obs[system_key] = system_name
    
    adata = sc.concat(list(adatas.values()))
    
    # Find common genes
    for name, ad in adatas.items():
        if name == list(adatas.keys())[0]:
            common_genes = set(ad.var_names)
        else:
            common_genes = common_genes.intersection(ad.var_names)
    
    adata = adata[:, list(common_genes)].copy()
    print(f"Common genes: {len(common_genes)}")
    
    # Store counts
    adata.layers["counts"] = adata.X.copy()
    
    # HVG selection
    sc.pp.highly_variable_genes(
        adata,
        n_top_genes=n_top_genes,
        flavor="seurat_v3",
        batch_key=system_key,
        layer="counts"
    )
    adata = adata[:, adata.var["highly_variable"]].copy()
    
    # Setup with system as covariate
    scvi.model.SCVI.setup_anndata(
        adata,
        layer="counts",
        batch_key=batch_key if batch_key in adata.obs else None,
        categorical_covariate_keys=[system_key]
    )
    
    # Train
    model = scvi.model.SCVI(adata, n_latent=n_latent, n_layers=2)
    model.train(max_epochs=300, early_stopping=True)
    
    # Get representation
    adata.obsm["X_integrated"] = model.get_latent_representation()
    
    # Clustering
    sc.pp.neighbors(adata, use_rep="X_integrated")
    sc.tl.umap(adata)
    sc.tl.leiden(adata)
    
    return adata, model

# Usage
adatas = {
    "10x_v3": sc.read_h5ad("10x_v3_data.h5ad"),
    "Smart-seq2": sc.read_h5ad("smartseq_data.h5ad"),
    "Drop-seq": sc.read_h5ad("dropseq_data.h5ad")
}

adata_integrated, model = integrate_cross_system(adatas)

# Visualize
sc.pl.umap(adata_integrated, color=["system", "leiden"])

Troubleshooting

Key References

• Lopez et al. (2018) "Deep generative modeling for single-cell transcriptomics"
• Luecken et al. (2022) "Benchmarking atlas-level data integration in single-cell genomics"

CITE-seq Analysis with totalVI

This reference covers multi-modal analysis of CITE-seq data (RNA + surface proteins) using totalVI.

Overview

CITE-seq combines:

• scRNA-seq (transcriptome)
• Protein surface markers (antibody-derived tags, ADT)

totalVI jointly models both modalities to:

• Integrate across batches
• Denoise protein signal
• Learn joint latent representation
• Enable cross-modal imputation

Prerequisites

import scvi
import scanpy as sc
import mudata as md
import numpy as np
import pandas as pd

print(f"scvi-tools version: {scvi.__version__}")

Step 1: Load CITE-seq Data

From 10x Genomics (Cell Ranger)

# 10x outputs separate gene expression and feature barcoding
adata_rna = sc.read_10x_h5("filtered_feature_bc_matrix.h5", gex_only=False)

# Separate RNA and protein
adata_protein = adata_rna[:, adata_rna.var['feature_types'] == 'Antibody Capture'].copy()
adata_rna = adata_rna[:, adata_rna.var['feature_types'] == 'Gene Expression'].copy()

print(f"RNA: {adata_rna.shape}")
print(f"Protein: {adata_protein.shape}")

From MuData

# If data is in MuData format
mdata = md.read_h5mu("cite_seq.h5mu")

adata_rna = mdata['rna'].copy()
adata_protein = mdata['protein'].copy()

Combine into Single AnnData

# totalVI expects protein data in obsm
adata = adata_rna.copy()

# Add protein expression to obsm
adata.obsm["protein_expression"] = adata_protein.X.toarray() if hasattr(adata_protein.X, 'toarray') else adata_protein.X

# Store protein names
adata.uns["protein_names"] = list(adata_protein.var_names)

Step 2: Quality Control

RNA QC

# Standard RNA QC
# Handle both human (MT-) and mouse (mt-, Mt-) mitochondrial genes
    adata.var['mt'] = (
        adata.var_names.str.startswith('MT-') |
        adata.var_names.str.startswith('mt-') |
        adata.var_names.str.startswith('Mt-')
    )
sc.pp.calculate_qc_metrics(adata, qc_vars=['mt'], inplace=True)

# Filter cells
adata = adata[adata.obs['n_genes_by_counts'] > 200].copy()
adata = adata[adata.obs['pct_counts_mt'] < 20].copy()

# Filter genes
sc.pp.filter_genes(adata, min_cells=3)

Protein QC

# Protein QC
protein_counts = adata.obsm["protein_expression"]
print(f"Protein counts per cell: min={protein_counts.sum(1).min():.0f}, max={protein_counts.sum(1).max():.0f}")

# Check for isotype controls
# Isotype controls should have low counts
protein_names = adata.uns["protein_names"]
for i, name in enumerate(protein_names):
    if 'isotype' in name.lower() or 'control' in name.lower():
        print(f"{name}: mean={protein_counts[:, i].mean():.1f}")

Step 3: Data Preparation

Store Raw Counts

# Store RNA counts
adata.layers["counts"] = adata.X.copy()

# Protein must be raw ADT counts (NOT CLR-normalized)
# WARNING: If importing from Seurat, ensure you use raw counts, not CLR-normalized data
# Seurat's NormalizeData(normalization.method = "CLR") transforms counts - use the original assay

HVG Selection for RNA

# Select HVGs for RNA
# Note: totalVI uses all proteins regardless of HVG

sc.pp.highly_variable_genes(
    adata,
    n_top_genes=4000,  # Use more for CITE-seq
    flavor="seurat_v3",
    batch_key="batch" if "batch" in adata.obs else None,
    layer="counts"
)

# Subset to HVGs
adata = adata[:, adata.var["highly_variable"]].copy()

Step 4: Setup and Train totalVI

# Setup AnnData for totalVI
scvi.model.TOTALVI.setup_anndata(
    adata,
    layer="counts",
    protein_expression_obsm_key="protein_expression",
    batch_key="batch"  # Optional
)

# Create model
model = scvi.model.TOTALVI(
    adata,
    n_latent=20,
    latent_distribution="normal"  # or "ln" for log-normal
)

# Train
model.train(
    max_epochs=200,
    early_stopping=True,
    batch_size=128
)

# Check training
model.history['elbo_train'].plot()

Step 5: Get Latent Representation

# Joint latent space
adata.obsm["X_totalVI"] = model.get_latent_representation()

# Clustering and visualization
sc.pp.neighbors(adata, use_rep="X_totalVI")
sc.tl.umap(adata)
sc.tl.leiden(adata, resolution=1.0)

sc.pl.umap(adata, color=['leiden', 'batch'])

Step 6: Denoised Protein Expression

# Get denoised protein values
# This removes background noise from protein measurements

_, protein_denoised = model.get_normalized_expression(
    return_mean=True,
    transform_batch="batch1"  # Optional: normalize to specific batch
)

# Add to adata
adata.obsm["protein_denoised"] = protein_denoised

# Visualize denoised proteins
protein_names = adata.uns["protein_names"]
for i, protein in enumerate(protein_names[:5]):
    adata.obs[f"denoised_{protein}"] = protein_denoised[:, i]

sc.pl.umap(adata, color=[f"denoised_{p}" for p in protein_names[:5]])

Step 7: Normalized RNA Expression

# Get normalized RNA expression
rna_normalized, _ = model.get_normalized_expression(
    return_mean=True
)

# Store
adata.layers["totalVI_normalized"] = rna_normalized

Step 8: Differential Expression

RNA Differential Expression

# DE between clusters
de_rna = model.differential_expression(
    groupby="leiden",
    group1="0",
    group2="1"
)

# Filter significant genes
de_sig = de_rna[
    (de_rna['is_de_fdr_0.05']) &
    (abs(de_rna['lfc_mean']) > 1)
]

print(f"Significant DE genes: {len(de_sig)}")

Protein Differential Expression

# Protein DE
de_protein = model.differential_expression(
    groupby="leiden",
    group1="0",
    group2="1",
    mode="protein"
)

print(de_protein.head(20))

Step 9: Visualization

Protein Expression on UMAP

# Denoised protein on UMAP
import matplotlib.pyplot as plt

proteins_to_plot = ["CD3", "CD4", "CD8", "CD19", "CD14"]

fig, axes = plt.subplots(1, len(proteins_to_plot), figsize=(4*len(proteins_to_plot), 4))
for ax, protein in zip(axes, proteins_to_plot):
    idx = adata.uns["protein_names"].index(protein)
    sc.pl.umap(
        adata,
        color=adata.obsm["protein_denoised"][:, idx],
        ax=ax,
        title=protein,
        show=False
    )
plt.tight_layout()

Joint Heatmap

# Heatmap of top genes and proteins per cluster
sc.pl.dotplot(
    adata,
    var_names=de_sig.index[:20].tolist(),
    groupby="leiden",
    layer="totalVI_normalized"
)

Step 10: Cell Type Annotation

# Use both RNA and protein markers for annotation

# RNA markers
rna_markers = {
    'T cells': ['CD3D', 'CD3E'],
    'CD4 T': ['CD4'],
    'CD8 T': ['CD8A', 'CD8B'],
    'B cells': ['CD19', 'MS4A1'],
    'Monocytes': ['CD14', 'LYZ']
}

# Check denoised protein expression
for i, protein in enumerate(adata.uns["protein_names"]):
    if any(m in protein for m in ['CD3', 'CD4', 'CD8', 'CD19', 'CD14']):
        print(f"{protein}: cluster means")
        for cluster in adata.obs['leiden'].unique():
            mask = adata.obs['leiden'] == cluster
            mean_expr = adata.obsm["protein_denoised"][mask, i].mean()
            print(f"  Cluster {cluster}: {mean_expr:.2f}")

Complete Pipeline

def analyze_citeseq(
    adata_rna,
    adata_protein,
    batch_key=None,
    n_top_genes=4000,
    n_latent=20
):
    """
    Complete CITE-seq analysis with totalVI.
    
    Parameters
    ----------
    adata_rna : AnnData
        RNA expression (raw counts)
    adata_protein : AnnData
        Protein expression (raw counts)
    batch_key : str, optional
        Batch column in obs
    n_top_genes : int
        Number of HVGs
    n_latent : int
        Latent dimensions
        
    Returns
    -------
    Tuple of (processed AnnData, trained model)
    """
    import scvi
    import scanpy as sc
    
    # Ensure same cells
    common_cells = adata_rna.obs_names.intersection(adata_protein.obs_names)
    adata = adata_rna[common_cells].copy()
    adata_protein = adata_protein[common_cells].copy()
    
    # Add protein to obsm
    adata.obsm["protein_expression"] = adata_protein.X.toarray() if hasattr(adata_protein.X, 'toarray') else adata_protein.X
    adata.uns["protein_names"] = list(adata_protein.var_names)
    
    # RNA QC
    # Handle both human (MT-) and mouse (mt-, Mt-) mitochondrial genes
    adata.var['mt'] = (
        adata.var_names.str.startswith('MT-') |
        adata.var_names.str.startswith('mt-') |
        adata.var_names.str.startswith('Mt-')
    )
    sc.pp.calculate_qc_metrics(adata, qc_vars=['mt'], inplace=True)
    adata = adata[adata.obs['pct_counts_mt'] < 20].copy()
    sc.pp.filter_genes(adata, min_cells=3)
    
    # Store counts
    adata.layers["counts"] = adata.X.copy()
    
    # HVG selection
    sc.pp.highly_variable_genes(
        adata,
        n_top_genes=n_top_genes,
        flavor="seurat_v3",
        batch_key=batch_key,
        layer="counts"
    )
    adata = adata[:, adata.var["highly_variable"]].copy()
    
    # Setup totalVI
    scvi.model.TOTALVI.setup_anndata(
        adata,
        layer="counts",
        protein_expression_obsm_key="protein_expression",
        batch_key=batch_key
    )
    
    # Train
    model = scvi.model.TOTALVI(adata, n_latent=n_latent)
    model.train(max_epochs=200, early_stopping=True)
    
    # Get representations
    adata.obsm["X_totalVI"] = model.get_latent_representation()
    rna_norm, protein_denoised = model.get_normalized_expression(return_mean=True)
    adata.layers["totalVI_normalized"] = rna_norm
    adata.obsm["protein_denoised"] = protein_denoised
    
    # Clustering
    sc.pp.neighbors(adata, use_rep="X_totalVI")
    sc.tl.umap(adata)
    sc.tl.leiden(adata)
    
    return adata, model

# Usage
adata, model = analyze_citeseq(
    adata_rna,
    adata_protein,
    batch_key="batch"
)

# Visualize
sc.pl.umap(adata, color=['leiden', 'batch'])

Troubleshooting

Key References

• Gayoso et al. (2021) "Joint probabilistic modeling of single-cell multi-omic data with totalVI"

Data Preparation for scvi-tools

This reference covers how to properly prepare AnnData objects for use with scvi-tools models.

Overview

Proper data preparation is critical for scvi-tools. Key requirements:

1. Raw counts (not normalized)
2. Highly variable gene selection
3. Proper setup_anndata() call

Step 1: Load and Inspect Data

import scanpy as sc
import scvi
import numpy as np

# Load data
adata = sc.read_h5ad("data.h5ad")

# Check what's in adata.X
print(f"Shape: {adata.shape}")
print(f"X dtype: {adata.X.dtype}")
print(f"X contains integers: {np.allclose(adata.X.data, adata.X.data.astype(int))}")
print(f"X min: {adata.X.min()}, max: {adata.X.max()}")

Verify Raw Counts

# scvi-tools needs INTEGER counts
# If X appears normalized, check for raw counts

if hasattr(adata, 'raw') and adata.raw is not None:
    print("Found adata.raw")
    # Use raw counts
    adata = adata.raw.to_adata()
    
# Or check layers
if 'counts' in adata.layers:
    print("Found counts layer")
    # Will specify layer in setup_anndata

Step 2: Basic Filtering

# Filter cells (standard QC)
sc.pp.filter_cells(adata, min_genes=200)
sc.pp.filter_cells(adata, max_genes=5000)

# Calculate mito percent if not present
# Handle both human (MT-) and mouse (mt-, Mt-) mitochondrial genes
adata.var['mt'] = (
    adata.var_names.str.startswith('MT-') |
    adata.var_names.str.startswith('mt-') |
    adata.var_names.str.startswith('Mt-')
)
sc.pp.calculate_qc_metrics(adata, qc_vars=['mt'], inplace=True)
adata = adata[adata.obs['pct_counts_mt'] < 20].copy()

# Filter genes
sc.pp.filter_genes(adata, min_cells=3)

print(f"After filtering: {adata.shape}")

Step 3: Store Raw Counts

Critical: Always preserve raw counts before any normalization.

# Store raw counts in a layer
adata.layers["counts"] = adata.X.copy()

# Now you can normalize for other purposes (HVG selection)
# But scvi will use the counts layer

Step 4: Highly Variable Gene Selection

scvi-tools works best with 1,500-5,000 HVGs.

For Single-Batch Data

# Normalize for HVG selection only
adata_hvg = adata.copy()
sc.pp.normalize_total(adata_hvg, target_sum=1e4)
sc.pp.log1p(adata_hvg)

# Select HVGs
sc.pp.highly_variable_genes(
    adata_hvg,
    n_top_genes=2000,
    flavor="seurat"  # or "cell_ranger"
)

# Transfer HVG annotation
adata.var['highly_variable'] = adata_hvg.var['highly_variable']

For Multi-Batch Data (Recommended)

# Use seurat_v3 flavor with batch_key
# This selects genes variable across batches
sc.pp.highly_variable_genes(
    adata,
    n_top_genes=2000,
    flavor="seurat_v3",
    batch_key="batch",  # Your batch column
    layer="counts"      # Use raw counts
)

Subset to HVGs

# Subset to highly variable genes
adata = adata[:, adata.var['highly_variable']].copy()
print(f"After HVG selection: {adata.shape}")

Step 5: Setup AnnData

The setup_anndata() function registers data for the model.

Basic Setup

scvi.model.SCVI.setup_anndata(
    adata,
    layer="counts"  # Specify layer with raw counts
)

With Batch Information

scvi.model.SCVI.setup_anndata(
    adata,
    layer="counts",
    batch_key="batch"  # Column in adata.obs
)

With Cell Type Labels (for scANVI)

scvi.model.SCANVI.setup_anndata(
    adata,
    layer="counts",
    batch_key="batch",
    labels_key="cell_type"  # Column with cell type labels
)

With Continuous Covariates

scvi.model.SCVI.setup_anndata(
    adata,
    layer="counts",
    batch_key="batch",
    continuous_covariate_keys=["percent_mito", "n_genes"]
)

With Categorical Covariates

scvi.model.SCVI.setup_anndata(
    adata,
    layer="counts",
    batch_key="batch",
    categorical_covariate_keys=["donor", "technology"]
)

Multi-Modal Data Setup

CITE-seq (for totalVI)

# Protein data in adata.obsm
# RNA in adata.X, protein in separate matrix

# Add protein data
adata.obsm["protein_expression"] = protein_counts  # numpy array

# Setup for totalVI
scvi.model.TOTALVI.setup_anndata(
    adata,
    layer="counts",
    batch_key="batch",
    protein_expression_obsm_key="protein_expression"
)

Multiome RNA+ATAC (for MultiVI)

# RNA and ATAC in separate AnnData objects or MuData

import mudata as md

# If using MuData
mdata = md.read("multiome.h5mu")

scvi.model.MULTIVI.setup_mudata(
    mdata,
    rna_layer="counts",
    protein_layer=None,
    batch_key="batch",
    modalities={"rna": "rna", "accessibility": "atac"}
)

Complete Preparation Pipeline

For a complete preparation function, use prepare_adata() from scripts/model_utils.py:

from model_utils import prepare_adata

# Prepare data with QC, HVG selection, and layer setup
adata = prepare_adata(
    adata,
    batch_key="batch",
    n_top_genes=2000,
    min_genes=200,
    max_mito_pct=20
)

# Then setup for your model
import scvi
scvi.model.SCVI.setup_anndata(adata, layer="counts", batch_key="batch")

This function handles:

• Mitochondrial QC filtering
• Cell and gene filtering
• Storing counts in layer
• HVG selection (batch-aware if batch_key provided)
• Subsetting to HVGs

Checking Setup

# View registered data
print(adata.uns['_scvi_manager_uuid'])
print(adata.uns['_scvi_adata_minify_type'])

# For scVI
scvi.model.SCVI.view_anndata_setup(adata)

Common Issues and Solutions

Data Format Reference

Required

• adata.X or adata.layers["counts"]: Raw integer counts (sparse OK)
• adata.obs: Cell metadata DataFrame
• adata.var: Gene metadata DataFrame

Recommended

• adata.obs["batch"]: Batch/sample identifiers
• adata.var["highly_variable"]: HVG boolean mask

For scANVI

• adata.obs["labels"]: Cell type annotations
• Can include "Unknown" for unlabeled cells

Environment Setup for scvi-tools

This reference covers installation and environment configuration for scvi-tools.

Installation Options

Option 1: Conda Environment (Recommended)

# Create environment with GPU support
conda create -n scvi-env python=3.10
conda activate scvi-env

# Install scvi-tools
pip install scvi-tools

# For GPU acceleration (recommended for large datasets)
pip install torch --index-url https://download.pytorch.org/whl/cu118

# Common dependencies
pip install scanpy leidenalg

Option 2: Pip Only

# Create virtual environment
python -m venv scvi-env
source scvi-env/bin/activate  # Linux/Mac
# scvi-env\Scripts\activate   # Windows

# Install
pip install scvi-tools scanpy

Option 3: With Spatial Analysis Support

conda create -n scvi-spatial python=3.10
conda activate scvi-spatial

pip install scvi-tools scanpy squidpy

Option 4: With MuData Support (Multiome)

pip install scvi-tools mudata muon

Verify Installation

import scvi
import torch
import scanpy as sc

print(f"scvi-tools version: {scvi.__version__}")
print(f"scanpy version: {sc.__version__}")
print(f"PyTorch version: {torch.__version__}")
print(f"GPU available: {torch.cuda.is_available()}")

if torch.cuda.is_available():
    print(f"GPU device: {torch.cuda.get_device_name(0)}")
    print(f"GPU memory: {torch.cuda.get_device_properties(0).total_memory / 1e9:.1f} GB")

GPU Configuration

Check CUDA Version

nvidia-smi
nvcc --version

PyTorch CUDA Versions

Memory Management

import torch

# Clear GPU cache between models
torch.cuda.empty_cache()

# Monitor memory usage
print(f"Allocated: {torch.cuda.memory_allocated() / 1e9:.2f} GB")
print(f"Cached: {torch.cuda.memory_reserved() / 1e9:.2f} GB")

Common Issues

Jupyter Setup

# Install Jupyter kernel
pip install ipykernel
python -m ipykernel install --user --name scvi-env --display-name "scvi-tools"

# For interactive plots
pip install matplotlib seaborn

Recommended Package Versions

For reproducibility, pin versions:

pip install \
    scvi-tools>=1.0.0 \
    scanpy>=1.9.0 \
    anndata>=0.9.0 \
    torch>=2.0.0

Version Compatibility Guide

scvi-tools 1.x vs 0.x API Changes

The 1.x release introduced breaking changes. Key differences:

Migration from 0.x to 1.x

# OLD (0.x) - DEPRECATED
import scvi
scvi.data.setup_anndata(adata, layer="counts", batch_key="batch")
model = scvi.model.SCVI(adata)

# NEW (1.x) - CURRENT
import scvi
scvi.model.SCVI.setup_anndata(adata, layer="counts", batch_key="batch")
model = scvi.model.SCVI(adata)

Model-Specific Setup (1.x)

Each model has its own setup method:

# scVI
scvi.model.SCVI.setup_anndata(adata, layer="counts", batch_key="batch")

# scANVI
scvi.model.SCANVI.setup_anndata(adata, layer="counts", batch_key="batch", labels_key="cell_type")

# totalVI
scvi.model.TOTALVI.setup_anndata(adata, layer="counts", protein_expression_obsm_key="protein")

# MultiVI (uses MuData)
scvi.model.MULTIVI.setup_mudata(mdata, rna_layer="counts", atac_layer="counts")

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

# veloVI
scvi.external.VELOVI.setup_anndata(adata, spliced_layer="spliced", unspliced_layer="unspliced")

Minimum Version Requirements

Check Your Versions

import scvi
import scanpy as sc
import anndata
import torch

print(f"scvi-tools: {scvi.__version__}")
print(f"scanpy: {sc.__version__}")
print(f"anndata: {anndata.__version__}")
print(f"torch: {torch.__version__}")

# Check if using 1.x API
if hasattr(scvi.model.SCVI, 'setup_anndata'):
    print("Using scvi-tools 1.x API")
else:
    print("WARNING: Using deprecated 0.x API - please upgrade")

Known Compatibility Issues

Upgrading scvi-tools

# Upgrade to latest
pip install --upgrade scvi-tools

# Upgrade all dependencies
pip install --upgrade scvi-tools scanpy anndata torch

# If you have issues, clean install
pip uninstall scvi-tools
pip cache purge
pip install scvi-tools

Testing Installation

# Quick test with sample data
import scvi
import scanpy as sc

# Load test dataset
adata = scvi.data.heart_cell_atlas_subsampled()
print(f"Loaded test data: {adata.shape}")

# Setup and create model (quick test)
scvi.model.SCVI.setup_anndata(adata, layer="counts", batch_key="cell_source")
model = scvi.model.SCVI(adata, n_latent=10)
print("Model created successfully")

# Quick training test (1 epoch)
model.train(max_epochs=1)
print("Training works!")

Label Transfer and Reference Mapping with scANVI

This reference covers using scANVI for transferring cell type annotations from a reference atlas to query data.

Overview

Reference mapping (also called "label transfer") uses a pre-trained model on annotated reference data to predict cell types in new, unannotated query data. This is faster than re-clustering and more consistent across studies.

scANVI excels at this because it:

• Jointly embeds reference and query in shared space
• Transfers labels probabilistically
• Handles batch effects between reference and query

When to Use Reference Mapping

• Annotating new dataset using existing atlas
• Consistent annotation across multiple studies
• Speed: no need to re-cluster and manually annotate
• Quality: leverage expert-curated reference annotations

Workflow Options

1. Train new model: Train scANVI on reference, then map query
2. Use pre-trained model: Load existing model (e.g., from Model Hub)
3. scArches: Extend existing model with query data (preserves reference)

Option 1: Train scANVI on Reference

Step 1: Prepare Reference Data

import scvi
import scanpy as sc

# Load reference atlas
adata_ref = sc.read_h5ad("reference_atlas.h5ad")

# Check annotations
print(f"Reference cells: {adata_ref.n_obs}")
print(f"Cell types: {adata_ref.obs['cell_type'].nunique()}")
print(adata_ref.obs['cell_type'].value_counts())

# Ensure raw counts
adata_ref.layers["counts"] = adata_ref.raw.X.copy() if adata_ref.raw else adata_ref.X.copy()

# HVG selection
sc.pp.highly_variable_genes(
    adata_ref,
    n_top_genes=3000,
    flavor="seurat_v3",
    batch_key="batch" if "batch" in adata_ref.obs else None,
    layer="counts"
)
adata_ref = adata_ref[:, adata_ref.var["highly_variable"]].copy()

Step 2: Train scANVI on Reference

# First train scVI (unlabeled)
scvi.model.SCVI.setup_anndata(
    adata_ref,
    layer="counts",
    batch_key="batch"
)

scvi_ref = scvi.model.SCVI(adata_ref, n_latent=30)
scvi_ref.train(max_epochs=200)

# Initialize scANVI from scVI
scanvi_ref = scvi.model.SCANVI.from_scvi_model(
    scvi_ref,
    labels_key="cell_type",
    unlabeled_category="Unknown"
)

# Train scANVI
scanvi_ref.train(max_epochs=50)

# Save for later use
scanvi_ref.save("scanvi_reference_model/")

Step 3: Prepare Query Data

# Load query data
adata_query = sc.read_h5ad("query_data.h5ad")

# CRITICAL: Use same genes as reference
common_genes = adata_ref.var_names.intersection(adata_query.var_names)
print(f"Common genes: {len(common_genes)}")

# Subset query to reference genes
adata_query = adata_query[:, adata_ref.var_names].copy()

# Handle missing genes (set to 0)
missing_genes = set(adata_ref.var_names) - set(adata_query.var_names)
if missing_genes:
    # Add missing genes with zero expression
    import numpy as np
    from scipy.sparse import csr_matrix
    
    zero_matrix = csr_matrix((adata_query.n_obs, len(missing_genes)))
    # ... concat and reorder to match reference
    
# Store counts
adata_query.layers["counts"] = adata_query.X.copy()

Step 4: Map Query to Reference

# Prepare query data for mapping
scvi.model.SCANVI.prepare_query_anndata(adata_query, scanvi_ref)

# Create query model from reference
scanvi_query = scvi.model.SCANVI.load_query_data(
    adata_query,
    scanvi_ref
)

# Fine-tune on query (optional but recommended)
scanvi_query.train(
    max_epochs=100,
    plan_kwargs={"weight_decay": 0.0}
)

# Get predictions
adata_query.obs["predicted_cell_type"] = scanvi_query.predict()

# Get prediction probabilities
soft_predictions = scanvi_query.predict(soft=True)
adata_query.obs["prediction_score"] = soft_predictions.max(axis=1)

Step 5: Evaluate Predictions

# Confidence scores
print(f"Mean prediction confidence: {adata_query.obs['prediction_score'].mean():.3f}")

# Low confidence predictions
low_conf = adata_query.obs['prediction_score'] < 0.5
print(f"Low confidence cells: {low_conf.sum()} ({low_conf.mean()*100:.1f}%)")

# Visualize
sc.pp.neighbors(adata_query, use_rep="X_scANVI")
sc.tl.umap(adata_query)
sc.pl.umap(adata_query, color=['predicted_cell_type', 'prediction_score'])

Option 2: Use Pre-Trained Models

From Model Hub

# scvi-tools maintains models on HuggingFace
# Check: https://huggingface.co/scvi-tools

# Example: Load pre-trained model
from huggingface_hub import hf_hub_download

model_path = hf_hub_download(
    repo_id="scvi-tools/example-model",
    filename="model.pt"
)

# Load model
model = scvi.model.SCANVI.load(model_path, adata=adata_query)

From Published Atlas

# Many atlases provide pre-trained models
# Example workflow with CellTypist-style model

# Download reference model
# model = scvi.model.SCANVI.load("atlas_model/", adata=adata_query)

Option 3: scArches for Incremental Updates

scArches extends a reference model without retraining from scratch:

# Load existing reference model
scanvi_ref = scvi.model.SCANVI.load("reference_model/")

# Surgery: prepare for query integration
scanvi_ref.freeze_layers()

# Map query data
scvi.model.SCANVI.prepare_query_anndata(adata_query, scanvi_ref)
scanvi_query = scvi.model.SCANVI.load_query_data(adata_query, scanvi_ref)

# Train only query-specific parameters
scanvi_query.train(
    max_epochs=200,
    plan_kwargs={"weight_decay": 0.0}
)

Visualize Reference and Query Together

# Concatenate for joint visualization
adata_ref.obs["dataset"] = "reference"
adata_query.obs["dataset"] = "query"

# Get latent representations
adata_ref.obsm["X_scANVI"] = scanvi_ref.get_latent_representation()
adata_query.obsm["X_scANVI"] = scanvi_query.get_latent_representation()

# Combine
adata_combined = sc.concat([adata_ref, adata_query])

# Compute combined UMAP
sc.pp.neighbors(adata_combined, use_rep="X_scANVI")
sc.tl.umap(adata_combined)

# Plot
sc.pl.umap(
    adata_combined,
    color=["dataset", "cell_type", "predicted_cell_type"],
    ncols=2
)

Quality Control for Predictions

Confidence Filtering

# Filter predictions by confidence
confidence_threshold = 0.7

high_conf = adata_query[adata_query.obs['prediction_score'] >= confidence_threshold].copy()
low_conf = adata_query[adata_query.obs['prediction_score'] < confidence_threshold].copy()

print(f"High confidence: {len(high_conf)} ({len(high_conf)/len(adata_query)*100:.1f}%)")
print(f"Low confidence: {len(low_conf)} ({len(low_conf)/len(adata_query)*100:.1f}%)")

Marker Validation

# Validate predictions with known markers
markers = {
    'T cells': ['CD3D', 'CD3E'],
    'B cells': ['CD19', 'MS4A1'],
    'Monocytes': ['CD14', 'LYZ']
}

for ct, genes in markers.items():
    ct_cells = adata_query[adata_query.obs['predicted_cell_type'] == ct]
    if len(ct_cells) > 0:
        for gene in genes:
            if gene in adata_query.var_names:
                expr = ct_cells[:, gene].X.mean()
                print(f"{ct} - {gene}: {expr:.3f}")

Complete Pipeline

def transfer_labels(
    adata_ref,
    adata_query,
    cell_type_key="cell_type",
    batch_key=None,
    n_top_genes=3000,
    confidence_threshold=0.5
):
    """
    Transfer cell type labels from reference to query.
    
    Parameters
    ----------
    adata_ref : AnnData
        Annotated reference data
    adata_query : AnnData
        Unannotated query data
    cell_type_key : str
        Column with cell type annotations in reference
    batch_key : str, optional
        Batch column
    n_top_genes : int
        Number of HVGs
    confidence_threshold : float
        Minimum confidence for predictions
        
    Returns
    -------
    AnnData with predictions
    """
    import scvi
    import scanpy as sc
    
    # Prepare reference
    adata_ref = adata_ref.copy()
    adata_ref.layers["counts"] = adata_ref.X.copy()
    
    sc.pp.highly_variable_genes(
        adata_ref,
        n_top_genes=n_top_genes,
        flavor="seurat_v3",
        batch_key=batch_key,
        layer="counts"
    )
    adata_ref = adata_ref[:, adata_ref.var["highly_variable"]].copy()
    
    # Train reference model
    scvi.model.SCVI.setup_anndata(adata_ref, layer="counts", batch_key=batch_key)
    scvi_ref = scvi.model.SCVI(adata_ref, n_latent=30)
    scvi_ref.train(max_epochs=200)
    
    scanvi_ref = scvi.model.SCANVI.from_scvi_model(
        scvi_ref,
        labels_key=cell_type_key,
        unlabeled_category="Unknown"
    )
    scanvi_ref.train(max_epochs=50)
    
    # Prepare query
    adata_query = adata_query[:, adata_ref.var_names].copy()
    adata_query.layers["counts"] = adata_query.X.copy()
    
    # Map query
    scvi.model.SCANVI.prepare_query_anndata(adata_query, scanvi_ref)
    scanvi_query = scvi.model.SCANVI.load_query_data(adata_query, scanvi_ref)
    scanvi_query.train(max_epochs=100, plan_kwargs={"weight_decay": 0.0})
    
    # Get predictions
    adata_query.obs["predicted_cell_type"] = scanvi_query.predict()
    soft = scanvi_query.predict(soft=True)
    adata_query.obs["prediction_score"] = soft.max(axis=1)
    
    # Mark low confidence
    adata_query.obs["confident_prediction"] = adata_query.obs["prediction_score"] >= confidence_threshold
    
    # Add latent representation
    adata_query.obsm["X_scANVI"] = scanvi_query.get_latent_representation()
    
    return adata_query, scanvi_ref, scanvi_query

# Usage
adata_annotated, ref_model, query_model = transfer_labels(
    adata_ref,
    adata_query,
    cell_type_key="cell_type"
)

# Visualize
sc.pp.neighbors(adata_annotated, use_rep="X_scANVI")
sc.tl.umap(adata_annotated)
sc.pl.umap(adata_annotated, color=['predicted_cell_type', 'prediction_score'])

Troubleshooting

Key References

• Xu et al. (2021) "Probabilistic harmonization and annotation of single-cell transcriptomics data with deep generative models"
• Lotfollahi et al. (2022) "Mapping single-cell data to reference atlases by transfer learning"

Multiome Analysis with MultiVI

This reference covers joint RNA and ATAC-seq analysis from multiome experiments using MultiVI.

Overview

MultiVI is a deep generative model for analyzing multiome data (simultaneous RNA-seq and ATAC-seq from the same cells). It:

• Learns a joint latent representation across modalities
• Handles missing modalities (RNA-only or ATAC-only cells)
• Enables batch correction across experiments
• Supports imputation of missing modalities

Prerequisites

import scvi
import scanpy as sc
import mudata as md
import numpy as np

print(f"scvi-tools version: {scvi.__version__}")

Data Formats

Option 1: MuData (Recommended)

# Load multiome data as MuData
mdata = md.read("multiome.h5mu")

# Structure:
# mdata.mod['rna']  - AnnData with RNA counts
# mdata.mod['atac'] - AnnData with ATAC counts

print(f"RNA: {mdata.mod['rna'].shape}")
print(f"ATAC: {mdata.mod['atac'].shape}")

Option 2: Separate AnnData Objects

# Load separately
adata_rna = sc.read_h5ad("rna.h5ad")
adata_atac = sc.read_h5ad("atac.h5ad")

# Ensure same cells in same order
common_cells = adata_rna.obs_names.intersection(adata_atac.obs_names)
adata_rna = adata_rna[common_cells].copy()
adata_atac = adata_atac[common_cells].copy()

Step 1: Prepare RNA Data

# RNA preprocessing (standard scvi-tools pipeline)
adata_rna = mdata.mod['rna'].copy()

# Filter
sc.pp.filter_cells(adata_rna, min_genes=200)
sc.pp.filter_genes(adata_rna, min_cells=3)

# Store counts
adata_rna.layers["counts"] = adata_rna.X.copy()

# HVG selection
sc.pp.highly_variable_genes(
    adata_rna,
    n_top_genes=4000,
    flavor="seurat_v3",
    layer="counts",
    batch_key="batch"  # If multiple batches
)

# Subset to HVGs
adata_rna = adata_rna[:, adata_rna.var['highly_variable']].copy()

Step 2: Prepare ATAC Data

# ATAC preprocessing
adata_atac = mdata.mod['atac'].copy()

# Filter peaks
sc.pp.filter_genes(adata_atac, min_cells=10)

# Binarize accessibility
adata_atac.X = (adata_atac.X > 0).astype(np.float32)

# Select top accessible peaks (if too many)
if adata_atac.n_vars > 50000:
    peak_accessibility = np.array(adata_atac.X.sum(axis=0)).flatten()
    top_peaks = np.argsort(peak_accessibility)[-50000:]
    adata_atac = adata_atac[:, top_peaks].copy()

# Store in layer
adata_atac.layers["counts"] = adata_atac.X.copy()

Step 3: Create Combined MuData

# Ensure matching cells
common_cells = adata_rna.obs_names.intersection(adata_atac.obs_names)
adata_rna = adata_rna[common_cells].copy()
adata_atac = adata_atac[common_cells].copy()

# Create MuData
mdata = md.MuData({
    "rna": adata_rna,
    "atac": adata_atac
})

print(f"Combined multiome: {mdata.n_obs} cells")
print(f"RNA features: {mdata.mod['rna'].n_vars}")
print(f"ATAC features: {mdata.mod['atac'].n_vars}")

Step 4: Setup MultiVI

# Setup MuData for MultiVI
scvi.model.MULTIVI.setup_mudata(
    mdata,
    rna_layer="counts",
    atac_layer="counts",
    batch_key="batch",  # Optional
    modalities={
        "rna_layer": "rna",
        "batch_key": "rna",
        "atac_layer": "atac"
    }
)

Step 5: Train MultiVI

# Create model
model = scvi.model.MULTIVI(
    mdata,
    n_latent=20,
    n_layers_encoder=2,
    n_layers_decoder=2
)

# Train
model.train(
    max_epochs=300,
    early_stopping=True,
    early_stopping_patience=10,
    batch_size=128
)

# Check training
model.history['elbo_train'].plot()

Step 6: Get Joint Representation

# Latent representation
latent = model.get_latent_representation()

# Add to MuData
mdata.obsm["X_MultiVI"] = latent

# Clustering on joint space
sc.pp.neighbors(mdata, use_rep="X_MultiVI")
sc.tl.umap(mdata)
sc.tl.leiden(mdata, resolution=1.0)

# Visualize
sc.pl.umap(mdata, color=['leiden', 'batch'], ncols=2)

Step 7: Modality-Specific Analysis

Impute Missing Modality

# Impute RNA expression for ATAC-only cells
# (Useful when integrating with ATAC-only datasets)
imputed_rna = model.get_normalized_expression(
    modality="rna"
)

# Impute accessibility for RNA-only cells
imputed_atac = model.get_accessibility_estimates()

Differential Analysis

# Differential expression (RNA)
de_results = model.differential_expression(
    groupby="leiden",
    group1="0",
    group2="1"
)

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

Handling Partial Data

MultiVI can integrate datasets with only one modality:

# Dataset 1: Full multiome
# Dataset 2: RNA only
# Dataset 3: ATAC only

# Mark missing modalities
mdata.obs['modality'] = 'paired'  # For cells with both
# For RNA-only cells, ATAC data should be missing/NaN
# For ATAC-only cells, RNA data should be missing/NaN

# MultiVI handles this automatically during training

Complete Pipeline

def analyze_multiome(
    adata_rna,
    adata_atac,
    batch_key=None,
    n_top_genes=4000,
    n_top_peaks=50000,
    n_latent=20,
    max_epochs=300
):
    """
    Complete multiome analysis with MultiVI.

    Parameters
    ----------
    adata_rna : AnnData
        RNA count data
    adata_atac : AnnData
        ATAC peak data
    batch_key : str, optional
        Batch column name
    n_top_genes : int
        Number of HVGs for RNA
    n_top_peaks : int
        Number of top peaks for ATAC
    n_latent : int
        Latent dimensions
    max_epochs : int
        Maximum training epochs

    Returns
    -------
    MuData with joint representation
    """
    import scvi
    import scanpy as sc
    import mudata as md
    import numpy as np

    # Get common cells
    common_cells = adata_rna.obs_names.intersection(adata_atac.obs_names)
    adata_rna = adata_rna[common_cells].copy()
    adata_atac = adata_atac[common_cells].copy()

    # RNA preprocessing
    sc.pp.filter_genes(adata_rna, min_cells=3)
    adata_rna.layers["counts"] = adata_rna.X.copy()

    if batch_key:
        sc.pp.highly_variable_genes(
            adata_rna, n_top_genes=n_top_genes,
            flavor="seurat_v3", layer="counts", batch_key=batch_key
        )
    else:
        sc.pp.normalize_total(adata_rna, target_sum=1e4)
        sc.pp.log1p(adata_rna)
        sc.pp.highly_variable_genes(adata_rna, n_top_genes=n_top_genes)
        adata_rna.X = adata_rna.layers["counts"].copy()

    adata_rna = adata_rna[:, adata_rna.var['highly_variable']].copy()

    # ATAC preprocessing
    sc.pp.filter_genes(adata_atac, min_cells=10)
    adata_atac.X = (adata_atac.X > 0).astype(np.float32)

    if adata_atac.n_vars > n_top_peaks:
        peak_acc = np.array(adata_atac.X.sum(axis=0)).flatten()
        top_idx = np.argsort(peak_acc)[-n_top_peaks:]
        adata_atac = adata_atac[:, top_idx].copy()

    adata_atac.layers["counts"] = adata_atac.X.copy()

    # Create MuData
    mdata = md.MuData({"rna": adata_rna, "atac": adata_atac})

    # Setup and train
    scvi.model.MULTIVI.setup_mudata(
        mdata,
        rna_layer="counts",
        atac_layer="counts",
        batch_key=batch_key,
        modalities={"rna_layer": "rna", "batch_key": "rna", "atac_layer": "atac"}
    )

    model = scvi.model.MULTIVI(mdata, n_latent=n_latent)
    model.train(max_epochs=max_epochs, early_stopping=True)

    # Add representation
    mdata.obsm["X_MultiVI"] = model.get_latent_representation()

    # Cluster
    sc.pp.neighbors(mdata, use_rep="X_MultiVI")
    sc.tl.umap(mdata)
    sc.tl.leiden(mdata)

    return mdata, model


# Usage
mdata, model = analyze_multiome(
    adata_rna,
    adata_atac,
    batch_key="sample"
)

sc.pl.umap(mdata, color=['leiden', 'sample'])

Peak-to-Gene Linking

# Link ATAC peaks to genes based on correlation in latent space
# This identifies regulatory relationships

def link_peaks_to_genes(model, mdata, distance_threshold=100000):
    """
    Link peaks to nearby genes based on correlation.

    Parameters
    ----------
    model : MULTIVI
        Trained model
    mdata : MuData
        Multiome data
    distance_threshold : int
        Maximum distance (bp) to link peak to gene

    Returns
    -------
    DataFrame of peak-gene links
    """
    # Get imputed values
    rna_imputed = model.get_normalized_expression()
    atac_imputed = model.get_accessibility_estimates()

    # Correlate peak accessibility with gene expression
    # for peaks near gene promoters
    # ... (requires genomic coordinates)

    return peak_gene_links

Troubleshooting

Key References

• Ashuach et al. (2023) "MultiVI: deep generative model for the integration of multimodal data"

RNA Velocity with veloVI

This reference covers RNA velocity analysis using veloVI, a deep learning approach that improves upon traditional velocity methods.

Overview

RNA velocity estimates the future state of cells by modeling:

• Unspliced RNA: Newly transcribed, contains introns
• Spliced RNA: Mature mRNA, introns removed

The ratio of unspliced to spliced indicates whether a gene is being upregulated or downregulated.

Why veloVI?

Traditional methods (velocyto, scVelo) have limitations:

• Assume steady-state or dynamical model
• Sensitive to noise
• Don't handle batch effects

veloVI addresses these with:

• Probabilistic modeling
• Better uncertainty quantification
• Integration with scVI framework

Prerequisites

import scvi
import scvelo as scv
import scanpy as sc
import numpy as np

print(f"scvi-tools version: {scvi.__version__}")
print(f"scvelo version: {scv.__version__}")