Differential expression (DE) analysis¶

Compare protein abundances between G1 and G2 using the processed DDA-LFQ dataset, inspect the results, and create volcano plots.

For method selection and interpretation, see How it works: DE analysis.

Load required packages¶

In [1]:

import msmu as mm
import plotly.io as pio

pio.renderers.default = "png"

import msmu as mm import plotly.io as pio pio.renderers.default = "png"

Load processed data¶

Complete DDA-LFQ with the current msmu version, then run this notebook from the same working directory. Load its processed dda_lfq_PXD012986.h5mu result, which already contains filtered, log2-transformed, and normalized protein data. G1 and G2 each contain three samples.

In [2]:

mdata = mm.read_h5mu("dda_lfq_PXD012986.h5mu")

mdata = mm.read_h5mu("dda_lfq_PXD012986.h5mu")

Compare the two groups¶

Use mm.tl.run_de() to compare G2 (expr) against G1 (ctrl), using the factor value[condition] column in the protein modality's .obs. This call uses the default limma method and Benjamini–Hochberg correction (p_adjust="auto").

Set min_pct=0.5 to require observations in at least two of the three samples in each group. The input is already log2-transformed, so keep the default log_transformed=True; this argument describes the input and does not transform it. See the run_de() API reference for all parameters.

In [3]:

de_res = mm.tl.run_de(
    mdata,
    modality="protein",
    category="factor value[condition]",  # category in .obs to define groups
    ctrl="G1",                           # control group
    expr="G2",                           # experimental group
    min_pct=0.5,                         # require observations in at least 2 of 3 samples per group
)

de_res = mm.tl.run_de( mdata, modality="protein", category="factor value[condition]", # category in .obs to define groups ctrl="G1", # control group expr="G2", # experimental group min_pct=0.5, # require observations in at least 2 of 3 samples per group )

WARNING - mm.tl.run_de: the default DE engine is now 'limma' (previously a permutation test). Pass stat_method explicitly to choose an engine. This notice will be removed in a future release.

Inspect DE results¶

Use DeaResult.to_df() to obtain a pandas DataFrame. Preview the feature identifiers, fold changes, p-values, and q-values below. For this comparison, a positive log2fc means higher abundance in G2. Use q_value for significance decisions across the tested proteins; p_value is unadjusted.

pct_ctrl and pct_expr show the percentages of observed samples in each group. Features without a usable test result have NaN in p_value and q_value.

In [4]:

de_res.to_df()[["features", "log2fc", "p_value", "q_value"]].head()

de_res.to_df()[["features", "log2fc", "p_value", "q_value"]].head()

Out[4]:

	features	log2fc	p_value	q_value
0	A0A023T778,G3UZW7,P61327,Q9CQL1	0.108213	4.140496e-01	0.871968
1	A0A068CB13,A2QBC3,B2CNX8,B2CNX9,G3XSF9,I6VCW8,...	3.701067	1.104529e-08	0.000005
2	A0A087WNT1,A0A087WPE4,A0A087WQE6,P83940	0.006681	9.658937e-01	0.999956
3	A0A087WNT3,A0A087WNU9,A0A087WP64,A0A087WPE6,A0...	NaN	NaN	NaN
4	A0A087WNY6,A0A087WQA5,A0A087WQX8,A0A087WRP4,A0...	0.033887	8.063740e-01	0.975291

Visualize the comparison¶

The volcano plot shows log2fc against -log10(p_value), using the raw p-value cutoff 0.05 and the fc_pct_5 fold-change threshold by default. Read q_value from the result table for multiple-testing significance.

In [5]:

de_res.plot_volcano()

de_res.plot_volcano()

Optional: Run a permutation test¶

To use label permutation, set stat_method="welch" and request the number of label shuffles with n_resamples. This call uses empirical FDR correction by default; another correction can be selected with p_adjust.

Only 20 distinct three-versus-three label splits are available, so this call enumerates all of them. Increasing n_resamples cannot add new splits. The reported q-values depend on the data and correction method; treat the small-sample warning below as a design-level diagnostic rather than a universal q-value bound. See How it works: DE analysis for method details.

In [6]:

de_res_perm = mm.tl.run_de(
    mdata,
    modality="protein",
    category="factor value[condition]",
    ctrl="G1",
    expr="G2",
    stat_method="welch",    # "welch", "student", or "wilcoxon"
    n_resamples=1000,       # only 20 distinct splits exist here, so all are enumerated
)

de_res_perm = mm.tl.run_de( mdata, modality="protein", category="factor value[condition]", ctrl="G1", expr="G2", stat_method="welch", # "welch", "student", or "wilcoxon" n_resamples=1000, # only 20 distinct splits exist here, so all are enumerated )

WARNING - n_resamples=1000 exceeds the 20 distinct 3 vs 3 splits that exist, so every split is enumerated instead and the null is built from 20 permutations. Drawing more would only resample the same splits.

WARNING - 3 vs 3 with fdr='empirical' cannot produce q < 0.05: the null is built from only 20 permutations and the observed labelling is one of them, which floors q at 0.068 no matter how strong the effect. Use stat_method='limma' for this design.

Running Permutations:   0%|          | 0/20 [00:00<?, ?it/s]

Running Permutations: 100%|██████████| 20/20 [00:00<00:00, 2113.48it/s]