DDA - label-free¶

This tutorial demonstrates how to analyze label-free quantification (LFQ) data from data-dependent acquisition (DDA) using the msmu package.

Data preparation¶

The original dataset, described by Uszkoreit et al. (2022), is available as PXD012986. The search was performed with Sage v0.14.7 (Lazear (2023)).

For demonstration purposes, the example dataset was reduced to six samples and a total of 5,000 PSMs.

In [1]:

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"

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"

Load required packages¶

If you haven't installed the msmu package yet, please follow the installation guide.

In [2]:

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"

Read data¶

You can read data from various proteomics software outputs. Below are examples for Sage, MaxQuant, and FragPipe formats.

For this tutorial, we will use the Sage output as an example.

read_sage() function reads the Sage output files (lfq.tsv, results.sage.tsv) and creates modalities at MuData object.

In [3]:

# Sage format
mdata = mm.read_sage(identification_file=sage_idents, quantification_file=sage_quants, label="label_free")

# MaxQuant format
# mdata = mm.read_maxquant(identification_file="path/to/evidence.txt", label="label_free", acquisition="dda")

# FragPipe format
# mdata = mm.read_fragpipe(
#     identification_file="path/to/psm.tsv",
#     quantification_file="path/to/combined_modified_peptide.tsv",
#     label="label_free",
#     acquisition="dda",
# )

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

mdata.obs

# 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")

# Sage format mdata = mm.read_sage(identification_file=sage_idents, quantification_file=sage_quants, label="label_free") # MaxQuant format # mdata = mm.read_maxquant(identification_file="path/to/evidence.txt", label="label_free", acquisition="dda") # FragPipe format # mdata = mm.read_fragpipe( # identification_file="path/to/psm.tsv", # quantification_file="path/to/combined_modified_peptide.tsv", # label="label_free", # acquisition="dda", # ) # 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) mdata.obs # 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")

INFO - Reading SAGE Identification data: 1 file(s)

INFO - Reading SAGE Quantification data: 1 file(s)

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.

Out[3]:

	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]
QExHF04026	G1-1	Mus musculus	not available	C2C12	myoblast	CVCL_0188	C2C12	normal	1	QExHF04026	proteomic profiling by mass spectrometry	PXD012986	data-dependent acquisition	1	1	label free sample	Q Exactive HF	NT=Trypsin;AC=MS:1001251	G1
QExHF04028	G2-1	Mus musculus	not available	C2C12	myoblast	CVCL_0188	C2C12	normal	1	QExHF04028	proteomic profiling by mass spectrometry	PXD012986	data-dependent acquisition	1	1	label free sample	Q Exactive HF	NT=Trypsin;AC=MS:1001251	G2
QExHF04036	G1-2	Mus musculus	not available	C2C12	myoblast	CVCL_0188	C2C12	normal	2	QExHF04036	proteomic profiling by mass spectrometry	PXD012986	data-dependent acquisition	1	1	label free sample	Q Exactive HF	NT=Trypsin;AC=MS:1001251	G1
QExHF04038	G2-2	Mus musculus	not available	C2C12	myoblast	CVCL_0188	C2C12	normal	2	QExHF04038	proteomic profiling by mass spectrometry	PXD012986	data-dependent acquisition	1	1	label free sample	Q Exactive HF	NT=Trypsin;AC=MS:1001251	G2
QExHF04046	G1-3	Mus musculus	not available	C2C12	myoblast	CVCL_0188	C2C12	normal	3	QExHF04046	proteomic profiling by mass spectrometry	PXD012986	data-dependent acquisition	1	1	label free sample	Q Exactive HF	NT=Trypsin;AC=MS:1001251	G1
QExHF04048	G2-3	Mus musculus	not available	C2C12	myoblast	CVCL_0188	C2C12	normal	3	QExHF04048	proteomic profiling by mass spectrometry	PXD012986	data-dependent acquisition	1	1	label free sample	Q Exactive HF	NT=Trypsin;AC=MS:1001251	G2

Saving raw MuData object¶

After reading the search output, saving the MuData object as an h5mu file is recommended for future use.

In [4]:

mdata.write_h5mu("dda_lfq_PXD012986_raw.h5mu")

mdata.write_h5mu("dda_lfq_PXD012986_raw.h5mu")

Handling PSM level¶

Filtering - PSM¶

You can filter the data based on the column values, such as q-value. You can also filter on the flags the reader records, such as contaminant, or on string containment.

Filtering is split into two steps: first, you mark a filter condition using mm.pp.add_filter(), and then you apply the filter using mm.pp.apply_filter().

Here, we keep protein groups with q-value < 0.01.

In [5]:

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

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

INFO - Applying var filters for psm: ['q_value_lt_0.01']

Out[5]:

MuData object with n_obs × n_vars = 6 × 7914
  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:	'_log', 'sdrf'
  2 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', 'decoy', 'filter', 'decoy_filter'
      varm:	'search_result', 'filter'
      layers:	None
    peptide:	6 × 3578
      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:	'level'
      layers:	None

Handling peptide level¶

Summarization - peptide¶

You can summarize PSM-level data to peptide-level data using the mm.pp.to_peptide() function.

In [6]:

mdata = mm.pp.to_peptide(mdata)

mdata = mm.pp.to_peptide(mdata)

INFO - Peptide-level identifications: 3683 (3664 at 1% FDR)

INFO - Using existing peptide quantification data.

Optional: Quantification with FlashLFQ¶

Peptide intensities can also be obtained with FlashLFQ (Millikin et al. (2018)). Export the filtered PSMs using mm.io.write_flashlfq_input(), then run FlashLFQ with the corresponding spectra. Attach the resulting QuantifiedPeptides.tsv using mm.io.add_quant(). This replaces the peptide intensity matrix while retaining the identification annotations. Continue with peptide filtering, normalization, and protein inference below.

The default workflow uses Sage intensities; uncomment the cell below to use FlashLFQ.

1. Export the filtered PSMs.
2. Run FlashLFQ with the corresponding spectra. Follow the official FlashLFQ instructions for installation and settings. For this example, you can skip the external run and use the provided result for the same six samples.
3. Attach the peptide intensities from QuantifiedPeptides.tsv. Replace the example URL with your own result when running FlashLFQ on other data.

In [7]:

# Export the filtered PSMs for FlashLFQ.
# mm.io.write_flashlfq_input(mdata, "flashlfq_input.tsv")

# Run FlashLFQ externally, or use the provided result below.
# flashlfq_peptides = (
#     "https://raw.githubusercontent.com/bertis-informatics/msmu/refs/heads/main/data/flashlfq/QuantifiedPeptides.tsv"
# )

# Replace peptide intensities with the FlashLFQ result.
# mdata = mm.io.add_quant(mdata, quant_data=flashlfq_peptides, quant_tool="flashlfq")

# Export the filtered PSMs for FlashLFQ. # mm.io.write_flashlfq_input(mdata, "flashlfq_input.tsv") # Run FlashLFQ externally, or use the provided result below. # flashlfq_peptides = ( # "https://raw.githubusercontent.com/bertis-informatics/msmu/refs/heads/main/data/flashlfq/QuantifiedPeptides.tsv" # ) # Replace peptide intensities with the FlashLFQ result. # mdata = mm.io.add_quant(mdata, quant_data=flashlfq_peptides, quant_tool="flashlfq")

Filtering - peptide¶

In [8]:

mdata = mm.pp.add_filter(mdata, modality="peptide", column="q_value", keep="lt", value=0.01)
mdata = mm.pp.apply_filter(mdata, modality="peptide", on="var")

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", on="var") mdata

INFO - Applying var filters for peptide: ['q_value_lt_0.01']

Out[8]:

MuData object with n_obs × n_vars = 6 × 8000
  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:	'_log', 'sdrf'
  2 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', 'decoy', 'filter', 'decoy_filter'
      varm:	'search_result', 'filter'
      layers:	None
    peptide:	6 × 3664
      var:	'peptide', 'proteins', 'stripped_peptide', 'count_psm', 'PEP', 'q_value'
      uns:	'level', 'decoy', 'filter', 'decoy_filter'
      varm:	'filter'
      layers:	None

Normalization¶

Here, we log2 transform and normalize the data at the peptide level.

Median centering normalization is applied using mm.pp.normalize() function.

In [9]:

mdata = mm.pp.log2_transform(mdata, modality="peptide")
mdata = mm.pp.normalize(mdata, modality="peptide", method="median")

mdata = mm.pp.log2_transform(mdata, modality="peptide") mdata = mm.pp.normalize(mdata, modality="peptide", method="median")

Protein inference¶

You can infer protein-level data from peptide-level data using the mm.pp.infer_protein() function.

In [10]:

mdata = mm.pp.infer_protein(mdata)

mdata = mm.pp.infer_protein(mdata)

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

Handling protein level¶

Summarization - protein¶

You can summarize peptide-level data to protein-level data using the mm.pp.to_protein() function.

By default, all peptides within a protein group are used for protein group quantification aggregation. Set top_n (e.g., 3) to use only the top-ranked peptides.

In [11]:

mdata = mm.pp.to_protein(mdata, top_n=3, rank_method="median_intensity")

mdata = mm.pp.to_protein(mdata, top_n=3, rank_method="median_intensity")

INFO - Applying var filters for peptide: ['q_value_lt_0.01', 'peptide_type_eq_unique']

INFO - Protein-level identifications: 1501 (1475 at 1% FDR)

Filtering - protein¶

In [12]:

mdata = mm.pp.add_filter(mdata, modality="protein", column="q_value", keep="lt", value=0.01)
mdata = mm.pp.apply_filter(mdata, modality="protein", on="var")

mdata

mdata = mm.pp.add_filter(mdata, modality="protein", column="q_value", keep="lt", value=0.01) mdata = mm.pp.apply_filter(mdata, modality="protein", on="var") mdata

INFO - Applying var filters for protein: ['q_value_lt_0.01']

Out[12]:

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:	'_log', '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', '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', 'normalisation'
      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

Visualization¶

ID plot¶

In [13]:

mm.pl.plot_id(mdata, modality="protein", colorby="factor value[condition]").update_layout(width=800, height=400, title="Protein-level identifications by condition")

mm.pl.plot_id(mdata, modality="protein", colorby="factor value[condition]").update_layout(width=800, height=400, title="Protein-level identifications by condition")

Intensity distribution plot¶

In [14]:

mm.pl.plot_intensity(mdata, modality="protein", groupby="factor value[condition]").update_layout(width=800, height=400, title="Protein-level intensities by condition")

mm.pl.plot_intensity(mdata, modality="protein", groupby="factor value[condition]").update_layout(width=800, height=400, title="Protein-level intensities by condition")

Saving processed MuData object¶

You can save MuData object into an h5mu file. This allows you to easily reload the processed data in future sessions without repeating the entire analysis pipeline.

In [15]:

mdata.write_h5mu("dda_lfq_PXD012986.h5mu")

mdata.write_h5mu("dda_lfq_PXD012986.h5mu")

References¶

Uszkoreit, J., Barkovits, K., Pacharra, S., Pfeiffer, K., Steinbach, S., Marcus, K., and Eisenacher, M. Dataset containing physiological amounts of spike-in proteins into murine C2C12 background as a ground truth quantitative LC-MS/MS reference. Data in Brief 43, 108435 (2022).

Lazear, M. R. Sage: an open-source tool for fast proteomics searching and quantification at scale. Journal of Proteome Research 22(11), 3652–3659 (2023).

Millikin, R. J., Solntsev, S. K., Shortreed, M. R., and Smith, L. M. Ultrafast peptide label-free quantification with FlashLFQ. Journal of Proteome Research 17(1), 386–391 (2018). https://doi.org/10.1021/acs.jproteome.7b00608.