PyDESeq2 differential expression: ID mapping, DE testing, fold-change thresholding, and GSEA enrichment visualization in OmicVerse.
Use this skill when a user wants to reproduce the DESeq2 workflow showcased in 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.
import omicverse as ov and ov.style() to standardise visuals.ov.io.read(..., index_col=0, header=1)..bam from column names with [c.split('/')[-1].replace('.bam', '') for c in data.columns].ov.utils.download_geneid_annotation_pair()gene_id with gene symbols using ov.bulk.Matrix_ID_mapping(data, 'genesets/pair_<GENOME>.tsv').dds = ov.bulk.pyDEG(data) from the mapped counts.dds.drop_duplicates_index() and confirm success in logs.treatment_groups and control_groups lists that match column names exactly.dds.deg_analysis(treatment_groups, control_groups, method='DEseq2') to invoke PyDESeq2.dds.result.shape) and optionally filter low-expression genes, e.g. dds.result.loc[dds.result['log2(BaseMean)'] > 1].dds.foldchange_set(fc_threshold=-1, pval_threshold=0.05, logp_max=6) to auto-pick fold-change cutoffs.dds.plot_volcano(...) and summarise key genes.dds.plot_boxplot(genes=[...], treatment_groups=..., control_groups=..., figsize=(2, 3)).ov.utils.download_pathway_database() and load them through ov.utils.geneset_prepare.rnk = dds.ranking2gsea().gsea_obj = ov.bulk.pyGSEA(rnk, pathway_dict) and call gsea_obj.enrichment() to compute terms.gsea_obj.plot_enrichment(...) and GSEA curves with gsea_obj.plot_gsea(term_num=..., ...).# 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 treatment_groups + control_groups:
assert g in data.columns, f"Sample '{g}' not in count matrix columns: {list(data.columns)}"
gene_id mapping depends on species; direct them to download the correct genome pair when results look sparse.t_deseq2.ipynbsample/counts.txtreference.md