Bulk RNA-seq DESeq2 analysis with omicverse

Overview

Use this skill when a user wants to reproduce the DESeq2 workflow showcased in [tdeseq2.ipynb](../../omicverseguide/docs/Tutorials-bulk/t_deseq2.ipynb). It covers loading raw featureCounts matrices, mapping Ensembl IDs to symbols, running PyDESeq2 via ov.bulk.pyDEG, and exploring downstream enrichment plots.

Instructions

1. Import and format the expression matrix

- Call import omicverse as ov and ov.style() to standardise visuals. - Read tab-separated count data from featureCounts using ov.io.read(..., index_col=0, header=1). - Strip trailing .bam from column names with [c.split('/')[-1].replace('.bam', '') for c in data.columns].

1. Map gene identifiers

- Ensure the appropriate mapping pair exists by running ov.utils.downloadgeneidannotationpair(). - Replace geneid with gene symbols using ov.bulk.MatrixIDmapping(data, 'genesets/pair_<GENOME>.tsv').

1. Initialise the DEG object

- Create dds = ov.bulk.pyDEG(data) from the mapped counts. - Resolve duplicate gene names with dds.dropduplicatesindex() and confirm success in logs.

1. Define contrasts and run DESeq2

- Collect sample labels into treatmentgroups and controlgroups lists that match column names exactly. - Execute dds.deganalysis(treatmentgroups, control_groups, method='DEseq2') to invoke PyDESeq2.

1. Filter and tune thresholds

- Inspect result shape (dds.result.shape) and optionally filter low-expression genes, e.g. dds.result.loc[dds.result['log2(BaseMean)'] > 1]. - Set thresholds via dds.foldchangeset(fcthreshold=-1, pvalthreshold=0.05, logpmax=6) to auto-pick fold-change cutoffs.

1. Visualise differential genes

- Draw volcano plots with dds.plotvolcano(...) and summarise key genes. - Produce per-gene boxplots: dds.plotboxplot(genes=[...], treatmentgroups=..., controlgroups=..., figsize=(2, 3)).

1. Run enrichment analyses (optional)

- Download enrichment libraries using ov.utils.downloadpathwaydatabase() and load them through ov.utils.genesetprepare. - Rank genes for GSEA with rnk = dds.ranking2gsea(). - Instantiate gseaobj = ov.bulk.pyGSEA(rnk, pathwaydict) and call gseaobj.enrichment() to compute terms. - Plot enrichment bubble charts via gseaobj.plotenrichment(...) and GSEA curves with gseaobj.plotgsea(term_num=..., ...).

1. Defensive validation

``python # Before PyDESeq2: verify count matrix contains raw integers (not log-transformed) import numpy as np if hasattr(data, 'values'): sample = data.values.flatten()[:1000] else: sample = np.array(data).flatten()[:1000] if np.any(sample != sample.astype(int)): print("WARNING: Data may not be raw counts. PyDESeq2 requires integer counts, not log-transformed.") # Verify treatment/control groups match column names for g in treatmentgroups + controlgroups: assert g in data.columns, f"Sample '{g}' not in count matrix columns: {list(data.columns)}" ``

1. Troubleshooting

- If PyDESeq2 raises errors about size factors, remind users to provide raw counts (not log-transformed data). - gene_id mapping depends on species; direct them to download the correct genome pair when results look sparse. - Large pathway libraries may require raising recursion limits or filtering to the top N terms before plotting.

Examples

• "Run PyDESeq2 on treated vs control replicates and highlight the top enriched WikiPathways terms."
• "Filter DEGs to genes with log2(BaseMean) > 1, auto-select fold-change cutoffs, and create volcano and boxplots."
• "Generate the ranked gene list for GSEA and plot the enrichment curve for the top pathway."

References

• Tutorial notebook: [tdeseq2.ipynb](../../omicverseguide/docs/Tutorials-bulk/t_deseq2.ipynb)
• Sample featureCounts matrix: [sample/counts.txt](../../sample/counts.txt)
• Quick copy/paste commands: [reference.md](reference.md)