single-cell-multiomics

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