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

# 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

# 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

# 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

# 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

# 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

# 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.