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¶

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"

Data Preparation¶

Use the MuData object as mdata processed in tutorials DDA-TMT, DDA-LFQ, or DIA-LFQ. It should have a modality to test, in this case, protein.

In this tutorial, we will use the mdata processed in DDA-LFQ tutorial.

In [2]:

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", sdrf_file=sdrf
)

mdata = mm.pp.add_filter(mdata, modality="psm", column="q_value", keep="lt", value=0.01)
mdata = mm.pp.add_filter(mdata, modality="psm", column="proteins", keep="not_contains", value="contam_")
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

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", sdrf_file=sdrf ) mdata = mm.pp.add_filter(mdata, modality="psm", column="q_value", keep="lt", value=0.01) mdata = mm.pp.add_filter(mdata, modality="psm", column="proteins", keep="not_contains", value="contam_") 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 - Applying var filters for psm: ['q_value_lt_0.01', 'proteins_not_contains_contam_']
INFO - Peptide-level identifications: 3634 (3615 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: 3651
INFO - Removed indistinguishable: 1586
INFO - Removed subsettable: 539
INFO - Removed subsumable: 2
INFO - Total protein groups: 1524
INFO - Applying var filters for peptide: ['q_value_lt_0.01', 'peptide_type_eq_unique']
INFO - Protein-level identifications: 1489 (1463 at 1% FDR)
INFO - Applying var filters for protein: ['q_value_lt_0.01']

Out[2]:

MuData object with n_obs × n_vars = 6 × 9337
  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_data_file', 'comment_fraction_identifier', 'comment_technical_replicate', 'comment_label', 'comment_instrument', 'comment_cleavage_agent_details', 'factor_value_condition'
  uns:	'sdrf', '_cmd'
  3 modalities
    psm:	6 x 4259
      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_data_file', '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', 'sdrf', 'filter', 'decoy_filter'
      varm:	'search_result', 'filter'
    peptide:	6 x 3615
      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_data_file', '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'
    protein:	6 x 1463
      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_data_file', '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'

(Optional) Load MuData from local file¶

If MuData were saved locally, provide the path to the file in mm.read_h5mu() function.

In [3]:

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

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

Run DE analysis¶

Permutation-based DE analysis basically provided by 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 allows you to choose the statistical test method, such as welch, student, or wilcoxon.

measure defines the measure of central tendency used for calculating fold-change, either median or mean.

min_pct sets the minimum fraction of non-missing values required in at least one group for a feature to be considered in the analysis.

n_resamples sets the number of resampling iterations for permutation testing; if set to None, simple parametric p-values are computed.

fdr specifies the multiple testing correction method, which can be empirical, bh (Benjamini-Hochberg), or None.

log_transformed indicates whether the data is log-transformed.

More explanation for run_de() can be found in the DE Analysis documentation and the mm.tl.run_de() API reference.

In an example below, we assign 1000 resamples, but only 20 permutations are performed because each group has 3 samples and 3 by 3 gives 20 combinations, which is exact test.

In [4]:

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
    stat_method="welch",  # statistical test method: "welch", "student", "wilcoxon", default "welch"
    measure="median",  # measure of central tendancy for calculating fold-change: "median", "mean", default "median"
    min_pct=0.5,  # minimum fraction of non-missing values in at least one group, default 0.5
    fdr="empirical",  # multiple testing correction method: "empirical", "bh", or None, default "empirical"
    n_resamples=1000,  # by default, 1000 resamples; if None, simple parametric p-values are computed
    log_transformed=True,  # whether data is log-transformed, default True
)

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 stat_method="welch", # statistical test method: "welch", "student", "wilcoxon", default "welch" measure="median", # measure of central tendancy for calculating fold-change: "median", "mean", default "median" min_pct=0.5, # minimum fraction of non-missing values in at least one group, default 0.5 fdr="empirical", # multiple testing correction method: "empirical", "bh", or None, default "empirical" n_resamples=1000, # by default, 1000 resamples; if None, simple parametric p-values are computed log_transformed=True, # whether data is log-transformed, default True )

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

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.

In [5]:

de_res

de_res

Out[5]:

DeaResult(stat_method='welch', 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=(1463,), dtype=object), repr_ctrl=array([23.032452, 25.08355 , 25.420649, ..., 26.587545, 25.655592,
       22.321047], shape=(1463,), dtype=float32), repr_expr=array([23.090721, 28.73562 , 25.331665, ..., 26.365978, 25.997116,
       22.347776], shape=(1463,), dtype=float32), pct_ctrl=array([100., 100., 100., ..., 100., 100., 100.], shape=(1463,)), pct_expr=array([100., 100., 100., ..., 100., 100., 100.], shape=(1463,)), log2fc=array([ 0.0582695 ,  3.652071  , -0.08898354, ..., -0.22156715,
        0.34152412,  0.02672958], shape=(1463,), dtype=float32))

In [6]:

de_res.to_df().head()

de_res.to_df().head()

Out[6]:

	features	repr_ctrl	repr_expr	pct_ctrl	pct_expr	log2fc	p_value	q_value
0	A0A023T778,G3UZW7,P61327,Q9CQL1	23.032452	23.090721	100.0	100.0	0.058270	0.443105	0.813088
1	A0A068CB13,A2QBC3,B2CNX8,B2CNX9,G3XSF9,I6VCW8,...	25.083549	28.735620	100.0	100.0	3.652071	0.000174	0.075433
2	A0A087WNT1,A0A087WPE4,A0A087WQE6,P83940	25.420649	25.331665	100.0	100.0	-0.088984	0.970467	0.904035
3	A0A087WNT3,A0A087WNU9,A0A087WP64,A0A087WPE6,A0...	NaN	NaN	0.0	0.0	NaN	NaN	NaN
4	A0A087WNY6,A0A087WQA5,A0A087WQX8,A0A087WRP4,A0...	25.827114	25.794582	100.0	100.0	-0.032532	0.821479	0.880609

After permutation test, guidance for log2FC threshold at 5% (or 1%) on two-sided tails of null distribution from permutations is also provided.

In [7]:

print(f"Log2FC threshold at 5%: {de_res.fc_pct_5}")

print(f"Log2FC threshold at 5%: {de_res.fc_pct_5}")

Log2FC threshold at 5%: 0.33

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).

In [8]:

de_res.plot_volcano()

de_res.plot_volcano()

In [9]:

de_res.plot_volcano(label_top=3)

de_res.plot_volcano(label_top=3)