Differential Expression (DE) Analysis¶
This tutorial demonstrates how to perform differential expression analysis using the msmu package. Details about DE analysis methods can be found in the Differential Expression (DE) Analysis section.
Load Required Packages¶
import msmu as mm
import plotly.io as pio
pio.renderers.default = "png"
base_dir = "https://raw.githubusercontent.com/bertis-informatics/msmu/refs/heads/dev/data/sage_lfq"
sage_idents = f"{base_dir}/sage/results.sage.tsv"
sage_quants = f"{base_dir}/sage/lfq.tsv"
sdrf = f"{base_dir}/meta.sdrf.tsv"
mdata = mm.read_sage(
identification_file=sage_idents,
quantification_file=sage_quants,
label="label_free",
)
# Sample metadata from the SDRF: keep the table on the container, then project it onto obs
mdata = mm.pp.attach_sdrf(mdata, sdrf)
mdata = mm.pp.apply_sdrf_to_obs(mdata)
# Without an SDRF, add_meta() joins any table (DataFrame, csv, tsv, parquet) onto obs instead:
# mdata = mm.pp.add_meta(mdata, sdrf, format="sdrf", metadata_on="assay name")
mdata = mm.pp.add_filter(mdata, modality="psm", column="q_value", keep="lt", value=0.01)
mdata = mm.pp.apply_filter(mdata, modality="psm")
mdata = mm.pp.to_peptide(mdata)
mdata = mm.pp.add_filter(mdata, modality="peptide", column="q_value", keep="lt", value=0.01)
mdata = mm.pp.apply_filter(mdata, modality="peptide")
mdata = mm.pp.log2_transform(mdata, modality="peptide")
mdata = mm.pp.normalise(mdata, modality="peptide", method="median")
mdata = mm.pp.infer_protein(mdata)
mdata = mm.pp.to_protein(mdata, top_n=3, rank_method="median_intensity")
mdata = mm.pp.add_filter(mdata, modality="protein", column="q_value", keep="lt", value=0.01)
mdata = mm.pp.apply_filter(mdata, modality="protein")
mdata
INFO - Validating SDRF metadata for https://raw.githubusercontent.com/bertis-informatics/msmu/refs/heads/dev/data/sage_lfq/meta.sdrf.tsv. INFO - SDRF validation succeeded for https://raw.githubusercontent.com/bertis-informatics/msmu/refs/heads/dev/data/sage_lfq/meta.sdrf.tsv. INFO - Applying var filters for psm: ['q_value_lt_0.01'] INFO - Peptide-level identifications: 3683 (3664 at 1% FDR) INFO - Using existing peptide quantification data. INFO - Applying var filters for peptide: ['q_value_lt_0.01'] INFO - Starting protein inference INFO - Initial proteins: 3721 INFO - Removed indistinguishable: 1624 INFO - Removed subsettable: 559 INFO - Removed subsumable: 2 INFO - Total protein groups: 1536 INFO - Applying var filters for peptide: ['q_value_lt_0.01', 'peptide_type_eq_unique'] INFO - Protein-level identifications: 1501 (1475 at 1% FDR) INFO - Applying var filters for protein: ['q_value_lt_0.01']
MuData object with n_obs × n_vars = 6 × 9475
obs: 'source name', 'characteristics[organism]', 'characteristics[organism part]', 'characteristics[cell line]', 'characteristics[cell type]', 'characteristics[cellosaurus accession]', 'characteristics[cellosaurus name]', 'characteristics[disease]', 'characteristics[biological replicate]', 'assay name', 'technology type', 'comment[proteomexchange accession number]', 'comment[proteomics data acquisition method]', 'comment[fraction identifier]', 'comment[technical replicate]', 'comment[label]', 'comment[instrument]', 'comment[cleavage agent details]', 'factor value[condition]'
uns: '_cmd', 'sdrf'
3 modalities
psm: 6 × 4336
obs: 'source name', 'characteristics[organism]', 'characteristics[organism part]', 'characteristics[cell line]', 'characteristics[cell type]', 'characteristics[cellosaurus accession]', 'characteristics[cellosaurus name]', 'characteristics[disease]', 'characteristics[biological replicate]', 'assay name', 'technology type', 'comment[proteomexchange accession number]', 'comment[proteomics data acquisition method]', 'comment[fraction identifier]', 'comment[technical replicate]', 'comment[label]', 'comment[instrument]', 'comment[cleavage agent details]', 'factor value[condition]'
var: 'proteins', 'peptide', 'stripped_peptide', 'filename', 'scan_num', 'charge', 'peptide_length', 'expmass', 'calcmass', 'rt', 'missed_cleavages', 'semi_enzymatic', 'contaminant', 'PEP', 'score', 'q_value', 'decoy'
uns: 'level', 'search_engine', 'quantification', 'label', 'acquisition', 'identification_file', 'quantification_file', 'decoy', 'filter', 'decoy_filter'
varm: 'search_result', 'filter'
layers: None
peptide: 6 × 3664
obs: 'source name', 'characteristics[organism]', 'characteristics[organism part]', 'characteristics[cell line]', 'characteristics[cell type]', 'characteristics[cellosaurus accession]', 'characteristics[cellosaurus name]', 'characteristics[disease]', 'characteristics[biological replicate]', 'assay name', 'technology type', 'comment[proteomexchange accession number]', 'comment[proteomics data acquisition method]', 'comment[fraction identifier]', 'comment[technical replicate]', 'comment[label]', 'comment[instrument]', 'comment[cleavage agent details]', 'factor value[condition]'
var: 'peptide', 'proteins', 'stripped_peptide', 'count_psm', 'PEP', 'q_value', 'protein_group', 'peptide_type'
uns: 'level', 'decoy', 'filter', 'decoy_filter'
varm: 'filter'
layers: None
protein: 6 × 1475
obs: 'source name', 'characteristics[organism]', 'characteristics[organism part]', 'characteristics[cell line]', 'characteristics[cell type]', 'characteristics[cellosaurus accession]', 'characteristics[cellosaurus name]', 'characteristics[disease]', 'characteristics[biological replicate]', 'assay name', 'technology type', 'comment[proteomexchange accession number]', 'comment[proteomics data acquisition method]', 'comment[fraction identifier]', 'comment[technical replicate]', 'comment[label]', 'comment[instrument]', 'comment[cleavage agent details]', 'factor value[condition]'
var: 'count_psm', 'count_stripped_peptide', 'PEP', 'q_value'
uns: 'level', 'decoy', 'filter', 'decoy_filter'
varm: 'filter'
layers: None
(Optional) Load MuData from local file¶
If MuData were saved locally, provide the path to the file in mm.read_h5mu() function.
# mdata = mm.read_h5mu("dda_lfq_PXD012986.h5mu")
Run DE analysis¶
DE analysis is provided by the mm.tl.run_de() function. Here, we will compare two conditions, G1 and G2, in the condition column of mdata.obs.
modality specifies which modality to perform DE analysis on. category indicates the column name in mdata.obs that contains the group labels.
ctrl and expr define the control and experimental groups, respectively.
stat_method selects the engine and defaults to limma, an empirical-Bayes moderated-t test that keeps power at small sample sizes. We omit it below so that default runs; the opt-in label-permutation engines (welch, student, wilcoxon) are demonstrated at the end of this tutorial.
min_pct sets the minimum non-missing coverage a feature needs in every group (not in at least one), applied as a count, max(1, ceil(min_pct * n)) per group. It defaults to 0.0 — the estimability floor, where every feature whose contrast can be estimated is tested — and is an opt-in stringency knob you raise to also demand a minimum coverage per group. We pass 0.5 here (2 of 3) as an explicit example; on complete data like this it changes nothing, but on sparse data a fixed 0.5 drops many features, disproportionately the significant ones, so choose it from the coverage you observe (pct_ctrl/pct_expr) rather than from which features come out significant. Features below the threshold are not dropped: they stay as rows with repr_*/pct_* filled but p_value/q_value set to NaN. On/off features (present in one group, absent in the other) are reported this way rather than tested — msmu does not impute — so read them off pct_ctrl/pct_expr.
log_transformed indicates whether the data is log-transformed.
The fold-change measure and the q-value method are decided by the engine, so there are no measure or fdr arguments: welch, student and limma are mean-based while wilcoxon is median-based, and the permutation engines report an empirical FDR while limma reports Benjamini-Hochberg q-values.
More explanation for run_de() can be found in the DE Analysis documentation and the mm.tl.run_de() API reference.
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, # (opt-in) also require >=50% coverage per group; default 0.0 tests every estimable feature
# log_transformed=True, # whether data is log-transformed, default True
) # stat_method omitted -> the default limma engine runs
WARNING:msmu._tools._dea: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.
Explore DE results¶
Result is stored in DeaResult object. And the object provides few methods to explore the results. This object provides to_df() method to convert the results into a pandas DataFrame for further exploration. And plot_volcano() method to visualize the results as a volcano plot.
de_res
DeaResult(stat_method='limma', ctrl='G1', expr='G2', features=array(['A0A023T778,G3UZW7,P61327,Q9CQL1',
'A0A068CB13,A2QBC3,B2CNX8,B2CNX9,G3XSF9,I6VCW8,P13006',
'A0A087WNT1,A0A087WPE4,A0A087WQE6,P83940', ..., 'Q9Z2U0', 'Q9Z2X1',
'Q9Z315'], shape=(1475,), dtype=object), repr_ctrl=array([23.02666752, 25.0174429 , 25.34271054, ..., 26.60428074,
25.71939277, 22.34628457], shape=(1475,)), repr_expr=array([23.1348809 , 28.71851014, 25.34939181, ..., 26.44222493,
26.00662334, 22.37461511], shape=(1475,)), pct_ctrl=array([100., 100., 100., ..., 100., 100., 100.], shape=(1475,)), pct_expr=array([100., 100., 100., ..., 100., 100., 100.], shape=(1475,)), log2fc=array([ 0.10821339, 3.70106724, 0.00668126, ..., -0.16205581,
0.28723057, 0.02833054], shape=(1475,)), contrast_label='G2 vs G1', fc_pct_1=0.62, fc_pct_5=0.27)
de_res.to_df().head()
| features | repr_ctrl | repr_expr | pct_ctrl | pct_expr | log2fc | statistic | p_value | q_value | |
|---|---|---|---|---|---|---|---|---|---|
| 0 | A0A023T778,G3UZW7,P61327,Q9CQL1 | 23.026668 | 23.134881 | 100.0 | 100.0 | 0.108213 | 0.864096 | 4.140496e-01 | 0.871968 |
| 1 | A0A068CB13,A2QBC3,B2CNX8,B2CNX9,G3XSF9,I6VCW8,... | 25.017443 | 28.718510 | 100.0 | 100.0 | 3.701067 | 25.988893 | 1.104529e-08 | 0.000005 |
| 2 | A0A087WNT1,A0A087WPE4,A0A087WQE6,P83940 | 25.342711 | 25.349392 | 100.0 | 100.0 | 0.006681 | 0.044192 | 9.658937e-01 | 0.999956 |
| 3 | A0A087WNT3,A0A087WNU9,A0A087WP64,A0A087WPE6,A0... | NaN | NaN | 0.0 | 0.0 | NaN | NaN | NaN | NaN |
| 4 | A0A087WNY6,A0A087WQA5,A0A087WQX8,A0A087WRP4,A0... | 25.826491 | 25.860378 | 100.0 | 100.0 | 0.033887 | 0.253835 | 8.063740e-01 | 0.975291 |
A fold-change guidance line is also provided: fc_pct_5 (and fc_pct_1) is the log2FC threshold at the 5% (or 1%) two-sided tail of a label-permutation null of the fold change. It is computed for every engine, so it is available here even though limma itself does not permute.
print(f"Log2FC threshold at 5%: {de_res.fc_pct_5}")
Log2FC threshold at 5%: 0.27
To visualize the DE results, plot_volcano() method is available. Log2FC and p-value thresholds can be set manually or automatically using the fc_pct_5 attribute from the result object and p-value of 0.05.
Top significant features can be labelled using label_top parameter (sorted by log2FC).
de_res.plot_volcano()
de_res.plot_volcano(label_top=3)
Permutation test (opt-in)¶
Passing welch, student, or wilcoxon to stat_method runs a non-parametric label-permutation test instead of limma. These engines always permute, and n_resamples sets the number of label shuffles.
Note the trade-off on a dataset this size. Each group has 3 samples, and 3 by 3 gives only 20 distinct label splits, so every split is enumerated (an exact test). Because the observed labelling is itself one of those 20, the empirical/BH q-value is floored at about 0.068 and cannot reach 0.05 no matter how strong the effect — run_de() emits a warning to that effect. That degeneracy at small sample sizes is exactly why limma is the default engine.
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, # number of label shuffles; only 20 splits exist here, so an exact test is used
)
# the permutation q-value cannot cross 0.05 for this 3-vs-3 design
print(f"smallest q-value from the permutation test: {de_res_perm.to_df().q_value.min():.3f}")
WARNING:msmu._statistics._permutation: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:msmu._statistics._permutation: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: 100%|██████████| 20/20 [00:00<00:00, 384.42it/s]
smallest q-value from the permutation test: 0.070