pydeseq2

import pandas as pd from pydeseq2.dds import DeseqDataSet from pydeseq2.ds import DeseqStats # 1. Create dummy count data (genes x samples) and metadata # In a real scenario, load with pd.read_csv('counts.csv', index_col=0).T # and pd.read_csv('metadata.csv', index_col=0) counts_df = pd.DataFrame( { 'sample1': [100, 50, 200, 10, 5, 150], 'sample2': [120, 60, 210, 12, 6, 160], 'sample3': [10, 5, 20, 100, 50, 15], 'sample4': [15, 7, 25, 110, 55, 18] }, index=['geneA', 'geneB', 'geneC', 'geneD', 'geneE', 'geneF'] ).T # Transpose to samples x genes metadata = pd.DataFrame( { 'condition': ['treated', 'treated', 'control', 'control'], 'batch': ['batch1', 'batch2', 'batch1', 'batch2'] }, index=['sample1', 'sample2', 'sample3', 'sample4'] ) # Ensure counts_df index (samples) matches metadata index (samples) assert counts_df.index.equals(metadata.index) # 2. Filter low-count genes (optional, but good practice) genes_to_keep = counts_df.columns[counts_df.sum(axis=0) >= 10] counts_df = counts_df[genes_to_keep] # 3. Initialize DeseqDataSet with a formulaic design dds = DeseqDataSet( counts=counts_df, metadata=metadata, design='~condition' ) # 4. Run the DESeq2 pipeline (normalization, dispersion, LFC estimation) dds.deseq2() # 5. Perform statistical testing deseq_stats = DeseqStats(dds, contrast=['condition', 'treated', 'control']) deseq_stats.wald_test() # 6. Apply LFC shrinkage (optional) deseq_stats.lfc_shrink(coeff='condition_treated_vs_control') # 7. Access results results = deseq_stats.results_df print(results.head()) # You can also access attributes directly from dds, e.g., normalized counts # print(dds.layers['normed_counts'].head())