Single Cell Multiomics

---
name: single-cell-multiomics
description: "Advanced single-cell multi-omics analysis including scRNA-seq, scCITE-seq, scATAC-seq, and TARGET-seq. Use when analyzing single-cell data, cell type identification, trajectory analysis, differential expression, UMAP/clustering, integrating protein and RNA modalities (TotalVI), or working with Scanpy, Seurat, scvi-tools. Includes workflows for MPN, hematologic malignancies, megakaryocyte biology."
license: Unknown
---

# Single-Cell Multi-Omics Analysis

## Choosing the Right scvi-tools Model (read FIRST)

Pick the model from the **measured modalities and the analysis goal** — do not
guess. The most common mistake is confusing TotalVI with MultiVI.

| Model | Use when… | Modalities | Notes |
|-------|-----------|------------|-------|
| **scVI** | Batch-correct / integrate / denoise scRNA-seq; get a latent space | RNA | Default for unimodal RNA integration |
| **scANVI** | scVI **and** you have (partial) cell-type labels to transfer/harmonize | RNA + labels | Semi-supervised label transfer |
| **TotalVI** | **CITE-seq**: paired RNA **and** surface protein (ADT) | RNA + protein | This is the RNA+protein model |
| **MultiVI** | **Multiome**: paired (or partially missing) ATAC **and** RNA | RNA + ATAC | Joint RNA+ATAC; imputes a missing modality |
| **PeakVI** | ATAC-only analysis (peak counts) | ATAC | Chromatin-accessibility latent space |
| **mrVI** | Compare many **samples/donors/conditions**; quantify sample-level effects | RNA (multi-sample) | Sample-stratified comparative analysis |

Decision tree:

1. Surface **protein (ADT)** measured alongside RNA (CITE-seq)? → **TotalVI**
2. **ATAC** measured alongside RNA (10x Multiome)? → **MultiVI** (use **PeakVI** for ATAC alone)
3. Comparing many samples/donors/conditions for sample-level effects? → **mrVI**
4. Have cell-type labels to transfer across batches? → **scANVI**
5. Otherwise — plain RNA batch integration / denoising? → **scVI**

> ⚠️ **Do not confuse these:** RNA + **protein** (CITE-seq) → **TotalVI**.
> RNA + **ATAC** (multiome / missing-modality imputation) → **MultiVI**.
> MultiVI is **not** for protein data. For multi-sample comparative analysis,
> reach for **mrVI**, not MultiVI/TotalVI.

## Core Libraries & Environment

```python
# Essential imports
import scanpy as sc
import anndata as ad
import scvi
import muon as mu
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns

# Settings
sc.settings.verbosity = 3
sc.settings.set_figure_params(dpi=100, frameon=False, figsize=(6, 6))
```

## Standard scRNA-seq Workflow

```python
# 1. Load and QC
adata = sc.read_10x_mtx('path/to/filtered_feature_bc_matrix/')
sc.pp.filter_cells(adata, min_genes=200)
sc.pp.filter_genes(adata, min_cells=3)
adata.var['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, :]

# 2. Normalization & HVG
sc.pp.normalize_total(adata, target_sum=1e4)
sc.pp.log1p(adata)
sc.pp.highly_variable_genes(adata, n_top_genes=2000, batch_key='batch')

# 3. Dimensionality reduction
sc.pp.scale(adata, max_value=10)
sc.tl.pca(adata, svd_solver='arpack')
sc.pp.neighbors(adata, n_neighbors=15, n_pcs=40)
sc.tl.umap(adata)
sc.tl.leiden(adata, resolution=0.5)
```

## TotalVI for CITE-seq Integration

```python
# Setup MuData
mdata = mu.MuData({'rna': adata_rna, 'protein': adata_prot})

# Train TotalVI
scvi.model.TOTALVI.setup_mudata(
    mdata, rna_layer='counts', protein_layer='counts',
    batch_key='batch', modalities={'rna_layer': 'rna', 'protein_layer': 'protein'}
)
model = scvi.model.TOTALVI(mdata, latent_distribution='normal', n_latent=20)
model.train(max_epochs=200, early_stopping=True)

# Get embeddings
mdata.obsm['X_totalVI'] = model.get_latent_representation()
sc.pp.neighbors(mdata, use_rep='X_totalVI')
sc.tl.umap(mdata)
sc.tl.leiden(mdata, key_added='leiden_totalVI', resolution=0.6)
```

## Differential Expression

```python
# DEG analysis
sc.tl.rank_genes_groups(adata, 'leiden', method='wilcoxon')
result = adata.uns['rank_genes_groups']
df = pd.DataFrame({
    'gene': result['names']['0'],
    'log2FC': result['logfoldchanges']['0'],
    'pval_adj': result['pvals_adj']['0']
})
sig_genes = df[(df['pval_adj'] < 0.05) & (abs(df['log2FC']) > 1)]
```

## Publication-Quality Visualization

```python
# Dot plot with proper expression cutoffs
sc.pl.dotplot(
    adata, var_names=marker_genes, groupby='leiden',
    expression_cutoff=0.0001, mean_only_expressed=False,
    standard_scale='None', smallest_dot=0.1, dot_max=1.0,
    cmap='viridis', colorbar_title='Expression'
)

# UMAP by batch
for batch in adata.obs['batch'].unique():
    adata_batch = adata[adata.obs['batch'] == batch]
    sc.pl.umap(adata_batch, color='FOXP3', title=f'{batch}')
```

## Cell Type Annotation Markers

### Hematopoietic Markers
- **HSC**: CD34, KIT, THY1, CD38low
- **CMP/GMP**: CD34+, CD38+, CD123
- **MEP**: CD34+, CD38+, CD41/ITGA2B
- **Megakaryocytes**: ITGA2B, PF4, GP1BA, PPBP, VWF
- **Erythroid**: HBB, HBA1/2, GYPA, KLF1

### MPN-Specific Markers
- **Inflammatory MKs**: S100A8/9, CHI3L1, CXCL8, IL6
- **Fibrosis markers**: TGFB1, COL1A1, LOXL2, VEGFA
- **Disease genes**: JAK2, CALR, MPL, PPM1D, ASXL1

## Output & Saving

```python
# Save processed data
adata.write('processed_adata.h5ad')
model.save('totalvi_model/')
df.to_csv('DEG_results.csv', index=False)
```

See `references/cell_markers.md` for complete marker lists.
See `references/scvi_advanced.md` for advanced scvi-tools workflows.