Bulk RNA-seq batch correction with ComBat

Overview

Apply this skill when a user has multiple bulk expression matrices measured across different batches and needs to harmonise them before downstream analysis. It follows [tbulkcombat.ipynb](../../omicverseguide/docs/Tutorials-bulk/tbulk_combat.ipynb), w hich demonstrates the pyComBat workflow on ovarian cancer microarray cohorts.

Instructions

1. Import core libraries

- Load omicverse as ov, anndata, pandas as pd, and matplotlib.pyplot as plt. - Call ov.ovplotset() (aliased ov.plot_set() in some releases) to align figures with omicverse styling.

1. Load each batch separately

- Read the prepared pickled matrices (or user-provided expression tables) with pd.readpickle(...)/pd.readcsv(...). - Transpose to gene × sample before wrapping them in anndata.AnnData objects so adata.obs stores sample metadata. - Assign a batch column for every cohort (adata.obs['batch'] = '1', '2', ...). Encourage descriptive labels when availa ble.

1. Concatenate on shared genes

- Use anndata.concat([adata1, adata2, adata3], merge='same') to retain the intersection of genes across batches. - Confirm the combined adata reports balanced sample counts per batch; if not, prompt users to re-check inputs.

1. Run ComBat batch correction

- Execute ov.bulk.batchcorrection(adata, batchkey='batch'). - Explain that corrected values are stored in adata.layers['batch_correction'] while the original counts remain in adata.X.

1. Export corrected and raw matrices

- Obtain DataFrames via adata.todf().T (raw) and adata.todf(layer='batchcorrection').T (corrected). - Encourage saving both tables (.tocsv(...)) plus the harmonised AnnData (adata.writeh5ad('adatabatch.h5ad', compressio n='gzip')).

1. Benchmark the correction

- For per-sample variance checks, draw before/after boxplots and recolour boxes using ov.pl.redcolor, bluecolor, gree ncolor palettes to match batches. - Copy raw counts to a named layer with adata.layers['raw'] = adata.X.copy() before PCA. - Run ov.pp.pca(adata, layer='raw', npcs=50) and ov.pp.pca(adata, layer='batchcorrection', npcs=50). - Visualise embeddings with ov.pl.embedding(..., basis='raw|original|X_pca', color='batch', frameon='small') and repeat fo r the corrected layer to verify mixing.

1. Defensive validation

``python # Before ComBat: verify batch column exists and has >1 batch assert 'batch' in adata.obs.columns, "adata.obs must contain a 'batch' column" nbatches = adata.obs['batch'].nunique() assert nbatches > 1, f"Only {nbatches} batch — need >1 for batch correction" # Verify gene overlap after concatenation if adata.nvars < 100: print(f"WARNING: Only {adata.n_vars} shared genes after concat — check gene ID harmonization") ``

1. Troubleshooting tips

- Mismatched gene identifiers cause dropped features—remind users to harmonise feature names (e.g., gene symbols) before conca tenation. - pyComBat expects log-scale intensities or similarly distributed counts; recommend log-transforming strongly skewed matrices. - If batchcorrection layer is missing, ensure the batchkey matches the column name in adata.obs.

Examples

• "Combine three GEO ovarian cohorts, run ComBat, and export both the raw and corrected CSV matrices."
• "Plot PCA embeddings before and after batch correction to confirm that batches 1–3 overlap."
• "Save the harmonised AnnData file so I can reload it later for downstream DEG analysis."

References

• Tutorial notebook: [tbulkcombat.ipynb](../../omicverseguide/docs/Tutorials-bulk/tbulk_combat.ipynb)
• Example inputs: [omicverseguide/docs/Tutorials-bulk/data/combat/](../../omicverseguide/docs/Tutorials-bulk/data/combat/)
• Quick copy/paste commands: [reference.md](reference.md)