Bulk RNA-seq differential expression with omicverse

Overview

Follow this skill to run the end-to-end differential expression (DEG) workflow showcased in [tdeg.ipynb](../../omicverseguide/docs/Tutorials-bulk/t_deg.ipynb). It assumes the user provides a raw gene-level count matrix (e.g., from featureCounts) and wants to analyse bulk RNA-seq cohorts inside omicverse.

Instructions

1. Set up the session

- Import omicverse as ov, scanpy as sc, and matplotlib.pyplot as plt. - Call ov.plot_set() so downstream plots adopt omicverse styling.

1. Prepare ID mapping assets

- When gene IDs must be converted to gene symbols, instruct the user to download mapping pairs via ov.utils.downloadgeneidannotation_pair() and store them under genesets/. - Mention the available prebuilt genomes (T2T-CHM13, GRCh38, GRCh37, GRCm39, danRer7, danRer11) and that users can generate their own mapping from GTF files if needed.

1. Load the raw counts

- Read tab-delimited featureCounts output with ov.pd.readcsv(..., sep='\t', header=1, indexcol=0). - Strip trailing .bam segments from column names using list comprehension so sample IDs are clean.

1. Map gene identifiers

- Run ov.bulk.MatrixIDmapping(countsdf, 'genesets/pair<GENOME>.tsv') to replace gene_id entries with gene symbols.

1. Initialise the DEG object

- Create dds = ov.bulk.pyDEG(mappedcounts). - Handle duplicate gene symbols with dds.dropduplicates_index() to keep the highest expressed version.

1. Normalise and estimate size factors

- Execute dds.normalize() to calculate DESeq2 size factors, correcting for library size and batch differences.

1. Run differential testing

- Collect treatment and control replicate labels into lists. - Call dds.deganalysis(treatmentgroups, control_groups, method='ttest') for the default Welch t-test. - Offer optional alternatives: method='edgepy' for edgeR-like tests and method='limma' for limma-style modelling.

1. Filter and threshold results

- Note that lowly expressed genes are retained by default; filter using dds.result.loc[dds.result['log2(BaseMean)'] > 1] when needed. - Set dynamic fold-change and significance cutoffs via dds.foldchangeset(fcthreshold=-1, pvalthreshold=0.05, logpmax=6) (fc_threshold=-1 auto-selects based on log2FC distribution).

1. Visualise differential expression

- Produce volcano plots with dds.plotvolcano(title=..., figsize=..., plotgenes=... or plotgenesnum=...) to highlight key genes. - Generate per-gene boxplots using dds.plotboxplot(genes=[...], treatmentgroups=..., controlgroups=..., figsize=..., legendbbox=...); adjust y-axis tick labels if required.

1. Perform pathway enrichment (optional)

- Download curated pathway libraries through ov.utils.downloadpathwaydatabase(). - Load genesets with ov.utils.genesetprepare(<path>, organism='Mouse'|'Human'|...). - Build the DEG gene list from dds.result.loc[dds.result['sig'] != 'normal'].index. - Run enrichment with ov.bulk.genesetenrichment(genelist=deggenes, pathwaysdict=..., pvaluetype='auto', organism=...). Encourage users without internet access to provide a background gene list. - Visualise single-library results via ov.bulk.genesetplot(...) and combine multiple ontologies using ov.bulk.genesetplotmulti(enrdict, colors_dict, num=...).

1. Document outputs

- Suggest exporting dds.result and enrichment tables to CSV for downstream reporting. - Encourage users to save figures generated by matplotlib (plt.savefig(...)) when running outside notebooks.

1. Defensive validation

``python # Before DEG: verify treatment/control groups exist as column names allcols = set(dds.result.columns) if hasattr(dds, 'result') else set(countsdf.columns) for g in treatmentgroups + controlgroups: assert g in allcols, f"Sample '{g}' not found in count matrix columns" # Verify groups don't overlap assert not set(treatmentgroups) & set(control_groups), "Treatment and control groups must not overlap" ``

1. Troubleshooting tips

- Ensure sample labels in treatmentgroups/controlgroups exactly match column names post-cleanup. - Verify required packages (omicverse, pyComplexHeatmap, gseapy) are installed for enrichment visualisations. - Remind users that internet access is required the first time they download gene mappings or pathway databases.

Examples

• "I have a featureCounts matrix for mouse tumour samples—normalize it with DESeq2, run t-test DEG, and highlight the top 8 genes in a volcano plot."
• "Use omicverse to compute edgeR-style differential expression between treated and control replicates, then run GO enrichment on significant genes."
• "Guide me through converting Ensembl IDs to symbols, performing limma DEG, and plotting boxplots for Krtap9-5 and Lef1."

References

• Detailed walkthrough notebook: [tdeg.ipynb](../../omicverseguide/docs/Tutorials-bulk/t_deg.ipynb)
• Sample count matrix for testing: [sample/counts.txt](../../sample/counts.txt)
• Quick copy/paste commands: [reference.md](reference.md)